c************************************************************************ c c SUBROUTINE DYICE c c Sea ice dynamic model: freedrift, cavitating, viscous plastic, granular c Steady state model : Acc term in the mom equ ignored c ==> No restriction on the time step c Advection scheme : upstream scheme c c Appropriate reference: c c Tremblay, L.-B. and Mysak, L.A., "Modelling sea ice as a granular c material including the dilatancy effect", JPO, (in press). c c Tremblay, L.-B. and Mysak, L.A., "The possible effects of including c ridge-related roughness in air-ice drag parameterization: a c sensitivity study", Annals of Gaciology, (in press). c c Flato, G.M. and Hibler III, W.D., 1992, "Modeling pack ice as a c cavitating fluid, JPO, Vol 22, 626-651. c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c************************************************************************ c======================================================================== c grid numbering is as follows: 'x' are extra grid cells used to c define the boundary conditions c c c B-grid C-grid c c c -|--------------|- -|--------------|- c | | | | c | | | | c variable location | Psi(i,j) | u(i,j) Psi(i,j) | c | | | | c u(i,j)| | | | c -|--------------|- -|----v(i,j)----|- c v(i,j) c c c |---|---|---|- |---|---|---|- c | x | | | | x | . | . | c |---|---|---|- |---|---|---|- c numbering | x | | | |0,1|1,1| . | c 0,1-1,1--|---|- |---|---|---|- c | x | x | x | |0,0|1,0| x | c 0,0-1,0--|---|- |---|---|---|- c c======================================================================== c************************************************************************* c subroutine dyice: c calculate the ice thickness (h), concentration (A) and ice velocity c (u_ice, v_ice) at the next time step c************************************************************************* subroutine dyice (nx,ny,npt,ntrac, + tracer_IO, ! tracers ( h = tracer(i,j,1) , A = tracer(i,j,2) ) + uice_IO, ! x-comp ice velocity + vice_IO, ! y-comp ice velocity + p_IO, ! internal ice pressure + etaC_IO, ! non-linear viscous coefficient + etaB_IO, ! non-linear viscous coefficient + uwnd, ! x-comp wind velocity + vwnd, ! y-comp wind velocity + u1, ! x-comp ocean surface current + v1, ! y-comp ocean surface current + tauwx, ! x-comp water drag + tauwy ! y-comp water drag + ) include 'comm_dyice.h' c nhn: rather than specifying A,B or C-grid, use {u,v,h,e}-points c where the h-points are at the cell centers with u,v arranged c as on a C-grid, and the e-points at the cell corners dimension + tracer_IO (npt,ntrac), ! passive tracer (h-point) + uice_IO (0:nx+2,0:ny+2),! x-comp ice velocity (u-point) + vice_IO (0:nx+2,0:ny+2),! y-comp ice velocity (v-point) + p_IO (0:nx+1,0:ny+1),! internal ice pressure (h-point) + etaC_IO (0:nx+1,0:ny+1),! non-linear viscous coeff (h-point) + etaB_IO (0:nx+2,0:ny+2),! non-linear viscous coeff (e-point) + uwnd (npt), ! x-comp wind velocity (h-point) + vwnd (npt), ! y-comp wind velocity (h-point) + u1 (npt), ! x-comp ocean current (h-point) + v1 (npt), ! y-comp ocean current (h-point) + tauwx (npt), ! x-comp water drag (h-point) + tauwy (npt) ! y-comp water drag (h-point) c------------------------------------------------------------------------ c non-dimensionalize the variables before solving for u_ice c------------------------------------------------------------------------ call convert ( nx,ny,ntrac,npt, + tracer_IO, uice_IO, vice_IO, p_IO, etaC_IO,etaB_IO, + uwnd, vwnd, u1, v1,tauwx, tauwy, 'dim2nondim' ) c------------------------------------------------------------------------ c solve for h and A (and other tracers) at the next time step c------------------------------------------------------------------------ call continuity (nx,ny,ntrac) c------------------------------------------------------------------------ c solve for the forcing terms (R) on the C-grid. c if InitialGuess = 1 (freedrift), it also provides a freedrift solution c if not, get an updated velocity field with the new forcing c------------------------------------------------------------------------ call forcing (nx,ny,ntrac) if ( InitialGuess .eq. 2) call UVsolve (nx,ny,ntrac) if ( InitialGuess .eq. 1) call freedrift (nx,ny,ntrac) call forcing (nx,ny,ntrac) c------------------------------------------------------------------------ c solve for the freedrift, cavitating or viscous plastic ice velocity c------------------------------------------------------------------------ if ( rheology .eq. 1) then do k = 1, 3 call freedrift (nx,ny,ntrac) call forcing (nx,ny,ntrac) enddo endif if ( rheology .eq. 2) then do k = 1, 3 call cavitating (nx,ny,ntrac) call forcing (nx,ny,ntrac) call UVsolve (nx,ny,ntrac) call forcing (nx,ny,ntrac) enddo endif if ( rheology .eq. 3) then do k = 1, 3 call cavitating (nx,ny,ntrac) call forcing (nx,ny,ntrac) call GranViscCoef (nx,ny,ntrac,k) call UVsolve (nx,ny,ntrac) call forcing (nx,ny,ntrac) enddo endif c------------------------------------------------------------------------ c dimensionalize the variables before leaving ice c------------------------------------------------------------------------ call convert ( nx,ny,ntrac,npt, + tracer_IO, uice_IO, vice_IO, p_IO, etaC_IO, etaB_IO, + uwnd, vwnd, u1, v1, tauwx, tauwy, 'nondim2dim' ) return end c************************************************************************ c Subroutine init_dyice: set constants and parameters c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c************************************************************************ subroutine init_dyice (nx,ny,ntrac,npt,iox,xdeg,ydeg,deltat_IO) include 'comm_dyice.h' parameter (REARTH = 6378000., * DTOSEC = 86400., * RAD = 3.14159265/180., * TOMEGA = 4.*3.14159265/DTOSEC) common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch dimension iox(npt),xdeg(nx),ydeg(ny) c Deltax grid size [m] c Deltay grid size [m] c Deltat time step [day] c ntrac number of ice tracers c Deltax = (blon-alon)*RAD/float(nx-1)*REARTH* cos(dyice_lat*rad) Deltax = (blon-alon)*RAD/float(nx-1)*REARTH !good for arctic only Deltay = (blat-alat)*RAD/float(ny-1)*REARTH dyice_fcor = TOMEGA * sin(dyice_lat*rad) Deltat = deltat_IO c------------------------------------------------------------------------ c set mask to 0 on land/bndy and to 1 in ocean using maskC c------------------------------------------------------------------------ nxy = nx * ny nxy1 = (nx+1)*(ny+1) nxy2 = (nx+2)*(ny+2) nxy3 = (nx+3)*(ny+3) call mem_alloc (p_h, nxy2, 2, 'h in dyice') call mem_alloc (p_A, nxy2, 2, 'A in dyice') call mem_alloc (p_tracer, nxy2*ntrac, 2, 'tracer in dyice') call mem_alloc (p_uice, nxy3, 2, 'uice in dyice') call mem_alloc (p_vice, nxy3, 2, 'vice in dyice') call mem_alloc (p_p, nxy2, 2, 'p in dyice') call mem_alloc (p_etaC, nxy2, 2, 'etaC in dyice') call mem_alloc (p_etaB, nxy3, 2, 'etaB in dyice') call mem_alloc (p_effshear, nxy2, 2, 'effshear in dyice') call mem_alloc (p_CdwA, nxy2, 2, 'CdwA in dyice') call mem_alloc (p_CdwC1, nxy1, 2, 'CdwC1 in dyice') call mem_alloc (p_CdwC2, nxy1, 2, 'CdwC2 in dyice') call mem_alloc (p_CdaA, nxy2, 2, 'CdaA in dyice') call mem_alloc (p_CdaC1, nxy1, 2, 'CdaC1 in dyice') call mem_alloc (p_CdaC2, nxy1, 2, 'CdaC2 in dyice') call mem_alloc (p_uair, nxy2, 2, 'uair in dyice') call mem_alloc (p_vair, nxy2, 2, 'vair in dyice') call mem_alloc (p_speeda, nxy2, 2, 'speeda in dyice') call mem_alloc (p_uwater, nxy2, 2, 'uwater in dyice') call mem_alloc (p_vwater, nxy2, 2, 'vwater in dyice') call mem_alloc (p_speediw, nxy2, 2, 'speediw in dyice') call mem_alloc (p_ufd, nxy2, 2, 'ufd in dyice') call mem_alloc (p_vfd, nxy2, 2, 'vfd in dyice') call mem_alloc1 (p_maskB, nxy3,'maskB in dyice') call mem_alloc1 (p_maskC, nxy2,'maskC in dyice') call mem_alloc1 (p_lsm, nxy,'lsm in dyice') call mem_alloc (p_R1, nxy3, 2, 'R1 in dyice') call mem_alloc (p_R2, nxy3, 2, 'R2 in dyice') call mem_alloc (p_x2, nxy2, 2, 'x2 in dyice') call mem_alloc (p_divflux, nxy, 2, 'divflux in dyice') call mem_alloc (p_R1p, nxy2, 2, 'R1p in dyice') call mem_alloc (p_R2p, nxy2, 2, 'R2p in dyice') call mem_alloc (p_R1pp, nxy2, 2, 'R1pp in dyice') call mem_alloc (p_R2pp, nxy2, 2, 'R2pp in dyice') call mem_alloc (p_A11, nxy2, 2, 'A11 in dyice') call mem_alloc (p_A12, nxy2, 2, 'A12 in dyice') call mem_alloc (p_A21, nxy2, 2, 'A21 in dyice') call mem_alloc (p_A22, nxy2, 2, 'A22 in dyice') call mem_alloc (p_detA, nxy2, 2, 'detA in dyice') call mem_alloc (p_uterm, nxy2, 2, 'uterm in dyice') call mem_alloc (p_vterm, nxy2, 2, 'vterm in dyice') call mem_alloc (p_dudy, nxy3, 2, 'dudy in dyice') call mem_alloc (p_dvdx, nxy3, 2, 'dvdx in dyice') call mem_alloc (p_Pmax, nxy1, 2, 'Pmax in dyice') call lsm_init (nx,ny,npt,iox,lsm) do j = 0, ny+1 do i = 0, nx+1 maskC(i,j) = 0 enddo enddo do j = 1, ny do i = 1, nx if (lsm(i,j).ne.0) maskC(i,j) = 1 enddo enddo do j = 1, ny if (maskC(1,j) .eq.1) maskC(0,j) = 1 if (maskC(nx,j).eq.1) maskC(nx+1,j) = 1 enddo do i = 1, nx if (maskC(i,1) .eq.1) maskC(i,0) = 1 if (maskC(i,ny).eq.1) maskC(i,ny+1) = 1 enddo do j = 0, ny+2 do i = 0, nx+2 maskB(i,j) = 0 enddo enddo do j = 1, ny+1 do i = 1, nx+1 maskB(i,j) = ( maskC(i,j) + maskC(i-1,j) + + maskC(i,j-1) + maskC(i-1,j-1) ) / 4 enddo enddo c------------------------------------------------------------------------ c Program constants c------------------------------------------------------------------------ pi = 4.0 * atan(1.) ! pi [] deg2rad = pi / 180.0 ! [rad / deg] day2sec = 24.0 * 60.0 ** 2 ! [s/day] c------------------------------------------------------------------------ c Material properties c------------------------------------------------------------------------ rhowater = 1.3 ! water density [kg/m3] rhoice = 9e02 ! ice density [kg/m3] rhoair = 1.3 ! air density [kg/m3] c------------------------------------------------------------------------ c Characteristic scales c------------------------------------------------------------------------ L = 1.6 ! char. length scale [m] Hi = 10.0 ! char. ice thickness scale [m] Ui = 1e-01 ! char. ice speed [m/s] Ua = 10.0 ! char. wind speed [m/s] c------------------------------------------------------------------------ c Characteristic quantities c------------------------------------------------------------------------ Taua = rhoair * Cdair * Ua ! char. wind stress (linear) Tauw = rhowater * Cdwater * Ui ! char. water drag (linear) if ( DragLaw .eq. 2) then Taua = Taua * Ua Tauw = Tauw * Ui elseif ( DragLaw .eq. 3) then Taua = Taua * Hi Tauw = Tauw * Hi elseif ( DragLaw .eq. 4) then Taua = Taua * Hi * Ua Tauw = Tauw * Hi * Ui endif Sigma = Pstar * Hi ! char. internal ice stress [N/m2] Fcor = rhoice* Hi* dyice_fcor* Ui ! char. Coriolis force [N/m2] TimeScale = L / Ui ! char. time [s] c------------------------------------------------------------------------ c Non-dimensional parameters c------------------------------------------------------------------------ etamax = etamax ! non-dim. max friction coeff + / (Sigma * L / Ui) do k = 1, ntrac Ktracer(k) = Ktracer(k) ! non-dim. diffusion coeff + / (Ui * L) enddo Deltat = Deltat * day2sec ! non-dim time step + / TimeScale Deltax = Deltax / L ! non-dim grid size Deltay = Deltay / L ! non-dim grid size rhof = Fcor / Taua ! non-dim Coriolis term Cdw = Tauw / Taua ! non-dim water drag gradp = Sigma / Taua / L phi = phi * deg2rad ! internal angle of friction [rad] delta = delta * deg2rad ! angle of dilatancy [rad] theta_a = theta_a * deg2rad ! wind turning angle [rad] theta_w = theta_w * deg2rad ! water turning angle [rad] sintheta_a = sin( theta_a ) costheta_a = cos( theta_a ) sintheta_w = sin( theta_w ) costheta_w = cos( theta_w ) tandelta = tan ( delta ) sinphi = sin (phi) return end c************************************************************************ c Subroutine convert: transfer from dimensional to nondimensional c variable and vice versa. c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c initial version c V02 19-05-97 N. Naik c revised for use with LOAM ocean model c c c************************************************************************ subroutine convert ( nx,ny,ntrac,npt, + tracer_IO, uice_IO, vice_IO, p_IO, etaC_IO, etaB_IO, + uwnd, vwnd, u1, v1, tauwx, tauwy, case ) include 'comm_dyice.h' parameter (TAUCON = 1030., TAUINV = 1./TAUCON) dimension + tracer_IO (npt,ntrac), ! passive tracer (h-points) + uice_IO (0:nx+2,0:ny+2), ! x-comp ice velocity (u-points) + vice_IO (0:nx+2,0:ny+2), ! y-comp ice velocity (v-points) + p_IO (0:nx+1,0:ny+1), ! internal ice pressure (h-points) + etaC_IO (0:nx+1,0:ny+1), ! non-linear viscous coeff (h-points) + etaB_IO (0:nx+2,0:ny+2), ! non-linear viscous coeff (e-points) + uwnd (npt), ! x-comp wind velocity (h-points) + vwnd (npt), ! y-comp wind velocity (h-points) + u1 (npt), ! x-comp ocean current (h-points) + v1 (npt), ! y-comp ocean current (h-points) + tauwx (npt), ! x-comp water drag (h-points) + tauwy (npt) ! y-comp water drag (h-points) character*10 case if ( case .eq. 'dim2nondim') then c------------------------------------------------------------------------ c non-dimensionalize the variables passed as arguments in 'ice' c------------------------------------------------------------------------ do j = 1, ny do i = 1, nx k = lsm(i,j) if (k.ne.0) then tracer(i,j,1) = tracer_IO(k,1)/ Hi do it = 2, ntrac tracer(i,j,it) = tracer_IO(k,it) enddo uair(i,j) = uwnd(k) / Ua vair(i,j) = vwnd(k) / Ua uwater(i,j) = u1(k) / Ui vwater(i,j) = v1(k) / Ui endif enddo enddo do j = 1, ny+1 do i = 1, nx+1 uice(i,j) = uice_IO(i,j) / Ui vice(i,j) = vice_IO(i,j) / Ui p(i,j) = p_IO(i,j) / Sigma etaC(i,j) = etaC_IO(i,j) * Ui / (Sigma * L) etaB(i,j) = etaB_IO(i,j) * Ui / (Sigma * L) enddo enddo elseif (case .eq. 'nondim2dim') then c------------------------------------------------------------------------ c dimensionalize the variables passed as arguments in 'ice' c------------------------------------------------------------------------ do j = 1, ny do i = 1, nx k = lsm(i,j) if (k.ne.0) then tracer_IO(k,1) = tracer(i,j,1) * Hi do it = 2, ntrac tracer_IO(k,it) = tracer(i,j,it) enddo twx= 0.5*(CdwC1(i,j)* uice(i,j)+ CdwC1(i+1,j)* uice(i+1,j)) twy= 0.5*(CdwC2(i,j)* vice(i,j)+ CdwC2(i,j+1)* vice(i,j+1)) tauwx(k) = TAUINV* twx * Tauw ! divide by density tauwy(k) = TAUINV* twy * Tauw endif uice_IO(i,j) = uice(i,j) * Ui vice_IO(i,j) = vice(i,j) * Ui p_IO (i,j) = p(i,j) * Sigma etaC_IO(i,j) = etaB(i,j) * ( Sigma * L ) / Ui etaB_IO(i,j) = etaB(i,j) * ( Sigma * L ) / Ui enddo enddo endif return end c*********************************************************************** c subroutine continuity: c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c*********************************************************************** subroutine continuity (nx,ny,ntrac) include 'comm_dyice.h' c------------------------------------------------------------------------ c maximum value for tracer #1 to 10 ( eg A <= 1 ). c------------------------------------------------------------------------ data tracermax /1e30, 1., 1e30, 1e30, 1e30, 1e30, 1e30, 1e30, + 1e30, 1e30/ data tracermin /0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0/ do k = 1, ntrac c------------------------------------------------------------------------ c compute the divergence of the flux of a given tracer c------------------------------------------------------------------------ do i = 1, nx do j = 1, ny dFxdx = ( fluxx(nx,ny,ntrac,tracer,i+1,j,k) - fluxx(nx,ny,ntrac,tracer,i,j,k) ) + / Deltax dFydy = ( fluxy(nx,ny,ntrac,tracer,i,j+1,k) - fluxy(nx,ny,ntrac,tracer,i,j,k) ) + / Deltay divflux(i,j) = dFxdx + dFydy enddo enddo c------------------------------------------------------------------------ c update the tracer values c (in a separate do-loop to conserve mass) c------------------------------------------------------------------------ do i = 1, nx do j = 1, ny tracer(i,j,k) = ( tracer(i,j,k) + - Deltat * divflux(i,j) ) * maskC(i,j) tracer(i,j,k) = min( max( tracer(i,j,k), tracermin(k) ), + tracermax(k) ) if ( k .eq. 1 ) h(i,j) = tracer(i,j,k) if ( k .eq. 2 ) A(i,j) = tracer(i,j,k) enddo enddo enddo c------------------------------------------------------------------------ c diffusion of a given tracer (if applicable) c------------------------------------------------------------------------ call diffusion ( nx,ny,ntrac,tracer ) return end c***************************************************************************** c Subroutine diffusion: c computes the diffusion of a given tracer from the (if applicable). c The same variable is used for the input and output. c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c***************************************************************************** subroutine diffusion (nx,ny,ntrac,x1) include 'comm_dyice.h' dimension x1(0:nx+1,0:ny+1,ntrac) integer tmpvar tiny = 1e-08 do k = 1, ntrac if ( Ktracer(k) .lt. tiny ) goto 30 factor = Ktracer(k) * Deltat / (Deltax ** 2) do i = 1, nx do j = 1, ny c------------------------------------------------------------------------ c d^2x/dx^2 c------------------------------------------------------------------------ tmpvar = maskC(i+1,j) - maskC(i-1,j) x2(i,j) = x1(i,j,k) if ( maskC(i,j) .eq. 0 ) goto 10 if (tmpvar .eq. 1) then x2(i,j) = x1(i,j,k) + + factor * (-x1(i,j,k) + x1(i+1,j,k) ) elseif (tmpvar .eq. -1) then x2(i,j) = x1(i,j,k) + + factor * (-x1(i,j,k) + x1(i-1,j,k) ) elseif (tmpvar .eq. 0) then x2(i,j) = x1(i,j,k) + + factor * ( x1(i-1,j,k) - 2.0 * x1(i,j,k) & + x1(i+1,j,k) ) endif 10 continue c------------------------------------------------------------------------ c d^2x/dy^2 c------------------------------------------------------------------------ tmpvar = maskC(i,j+1) - maskC(i,j-1) if ( maskC(i,j) .eq. 0 ) goto 20 if ( tmpvar .eq. 1 ) then x2(i,j) = x2(i,j) + + factor * (-x1(i,j,k) + x1(i,j+1,k) ) elseif ( tmpvar .eq.-1 ) then x2(i,j) = x2(i,j) + + factor * (-x1(i,j,k) + x1(i,j-1,k) ) elseif ( tmpvar .eq. 0 ) then x2(i,j) = x2(i,j) + + factor * ( x1(i,j-1,k) - 2.0 * x1(i,j,k) & + x1(i,j+1,k) ) endif 20 continue enddo enddo 30 continue enddo c------------------------------------------------------------------------ c set the output to the input variable c------------------------------------------------------------------------ do i = 1, nx do j = 1, ny x1(i,j,k) = x2(i,j) enddo enddo return end c**************************************************************************** c subroutine forcing: calculates the wind/ocean forcing c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c**************************************************************************** subroutine forcing (nx,ny,ntrac) include 'comm_dyice.h' do i = 0, nx+1 do j = 0, ny+1 k = lsm(i,j) if (k.ne.0) then c------------------------------------------------------------------------ c ice speed relative to ocean current and drag coefficients evaluated c at the grid center c------------------------------------------------------------------------ ucenter = ( uice(i,j) + uice(i+1,j) ) / 2.0 vcenter = ( vice(i,j) + vice(i,j+1) ) / 2.0 speediw(i,j) = sqrt( ( ucenter - uwater(i,j) ) ** 2 + + ( vcenter - vwater(i,j) ) ** 2 ) speeda(i,j) = sqrt( uair(i,j) ** 2 + vair(i,j) ** 2 ) if (DragLaw .eq. 1) then CdwA(i,j) = Cdw CdaA(i,j) = 1. * A(i,j) elseif (DragLaw .eq. 2) then CdwA(i,j) = Cdw * min(max(speediw(i,j), 6d-01), 3.0) CdaA(i,j) = 1. * A(i,j) * speeda(i,j) elseif (DragLaw .eq. 3) then CdwA(i,j) = Cdw * min(max(h(i,j), 0.1), 0.7) CdaA(i,j) = min(max(h(i,j), 0.1), 0.7) + * A(i,j) elseif (DragLaw .eq. 4) then CdwA(i,j) = Cdw * min(max(speediw(i,j), 6e-01), 3.0) + * min(max(h(i,j), 0.1), 0.7) CdaA(i,j) = 1. * speeda(i,j) + * min(max(h(i,j), 0.1), 0.7) + * A(i,j) endif endif enddo enddo c------------------------------------------------------------------------ c non-dimensional drag evaluated on the C-grid c------------------------------------------------------------------------ do i = 1, nx+1 do j = 1, ny CdwC1(i,j) = (CdwA(i,j) + CdwA(i-1,j)) / 2.0 enddo enddo do i = 1, nx do j = 1, ny+1 CdwC2(i,j) = (CdwA(i,j) + CdwA(i,j-1)) / 2.0 enddo enddo do i = 0, nx+1 do j = 0, ny+1 c------------------------------------------------------------------------ c Define the matrix [A] and vector {taua} component on the A-grid: c [A] {u_r} = {taua}, where u_r = u_i^fd - u_w^g. c------------------------------------------------------------------------ A11d = CdwA(i,j) * costheta_w A21d = rhof * h(i,j) + CdwA(i,j) * sintheta_w A12d = -A21d A22d = A11d taua1 = CdaA(i,j) * ( uair(i,j) * costheta_a + - vair(i,j) * sintheta_a ) taua2 = CdaA(i,j) * ( vair(i,j) * costheta_a + + uair(i,j) * sintheta_a ) c------------------------------------------------------------------------ c free drift ice velocities relative to geostr currents (u_r) (A-grid) c------------------------------------------------------------------------ detAd = A11d * A22d - A12d * A21d ufd(i,j) = ( A22d * taua1 - A12d * taua2 ) / detAd vfd(i,j) = (-A21d * taua1 + A11d * taua2 ) / detAd c------------------------------------------------------------------------ c free drift ice velocity (u_i^fd) (A-grid) c------------------------------------------------------------------------ ufd(i,j) = ufd(i,j) + uwater(i,j) vfd(i,j) = vfd(i,j) + vwater(i,j) enddo enddo do i = 1, nx+1 do j = 1, ny+1 c------------------------------------------------------------------------ c Define the wind stress, water drag and sea surface tilt on the C-grid: c {R} = [A] * {u_i^fd} c------------------------------------------------------------------------ A11d = CdwA(i,j) * costheta_w A21d = rhof * h(i,j) + CdwA(i,j) * sintheta_w A12d = -A21d A22d = A11d A2_11 = CdwA(i-1,j) * costheta_w A2_12 = -rhof * h(i-1,j) - CdwA(i-1,j) * sintheta_w A3_22 = CdwA(i,j-1) * costheta_w A3_21 = rhof * h(i,j-1) + CdwA(i,j-1) * sintheta_w R1(i,j) = ( A11d * ufd(i,j) + A12d * vfd(i,j) + + A2_11 * ufd(i-1,j) + A2_12 * vfd(i-1,j) ) / 2.0 R2(i,j) = ( A21d * ufd(i,j) + A22d * vfd(i,j) + + A3_21 * ufd(i,j-1) + A3_22 * vfd(i,j-1) ) / 2.0 enddo enddo return end c**************************************************************************** c subroutine UVsolve: c compute 'uice' and 'vice' at the next time step from the momentum c equations neglecting the acceleration term and keeping the new c pressure field and non-linear coefficient of friction constant. The c resulting set of equations ([A] {u_i,j} = R'') is linear and is c solved using Gauss-Seidel relaxation technique. c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c **************************************************************************** subroutine UVsolve (nx,ny,ntrac) include 'comm_dyice.h' integer uflag, vflag tol = 1e-03 ! tolerance on uice, vice nmax = 6000 ! max number of iterations k = 0 ! counter w = 0.7 ! relaxation coefficient c------------------------------------------------------------------------ c Forcing terms ( R' = R - grad(p) ) c------------------------------------------------------------------------ do i = 1, nx+1 do j = 1, ny+1 R1p(i,j) = R1(i,j) - ( p(i,j) - p(i-1,j) ) / Deltax R2p(i,j) = R2(i,j) - ( p(i,j) - p(i,j-1) ) / Deltay enddo enddo c------------------------------------------------------------------------ c Matrix [A]^{-1} components defined on the C-grid. c Calculated outside the loop since it is independent of 'uice'. c------------------------------------------------------------------------ do j = 1, ny+1 do i = 1, nx+1 A11tmp = etaB(i+1,j) + etaB(i,j) + etaC(i,j) + etaC(i,j-1) A11(i,j) = CdwC2(i,j) * costheta_w + A11tmp / Deltax ** 2 A22tmp = etaB(i,j+1) + etaB(i,j) + etaC(i,j) + etaC(i-1,j) A22(i,j) = CdwC1(i,j) * costheta_w + A22tmp / Deltax ** 2 uterm(i,j) = rhof * ( h(i,j) + h(i-1,j) ) / 2.0 + + CdwC1(i,j) * sintheta_w A12(i,j) = uterm(i,j) / 4.0 + ( etaC(i,j) - etaB(i,j) ) + / Deltax ** 2 vterm(i,j) = rhof * (h(i,j) + h(i,j-1)) / 2.0 + + CdwC2(i,j) * sintheta_w A21(i,j) =-vterm(i,j) / 4.0 + ( etaC(i,j) - etaB(i,j) ) + / Deltax ** 2 detA(i,j) = A11(i,j) * A22(i,j) - A12(i,j) * A21(i,j) enddo enddo c------------------------------------------------------------------------ c Solve for the ice velocity field, [A] {u_i,j} = R'' c------------------------------------------------------------------------ 5 continue error = 0.0 ! \ ierroru = 0 ! | Initialization of variables: jerroru = 0 ! > error between successive iterations & ierrorv = 0 ! | location (i,j) where the maximum error occurs jerrorv = 0 ! / do i = 1, nx+1 do j = 1, ny+1 upast = uice(i,j) vpast = vice(i,j) c------------------------------------------------------------------------ c forcing term R'' c------------------------------------------------------------------------ Sumv = ( vice(i,j+1) + vice(i-1,j) + vice(i-1,j+1) ) R1pptmp = ( etaC(i,j) * ( uice(i+1,j) - vice(i,j+1) ) + + etaC(i-1,j) * ( uice(i-1,j) + vice(i-1,j+1) + - vice(i-1,j) ) + + etaB(i,j+1) * ( uice(i,j+1) + vice(i,j+1) + - vice(i-1,j+1) ) + + etaB(i,j) * ( uice(i,j-1) + vice(i-1,j) ) ) + / Deltax ** 2 R1pp(i,j) = R1p(i,j) + uterm(i,j) * Sumv / 4.0 + R1pptmp Sumu = (uice(i+1,j) + uice(i,j-1) + uice(i+1,j-1)) R2pptmp = ( etaC(i,j) * ( vice(i,j+1) - uice(i+1,j) ) + + etaC(i,j-1) * ( vice(i,j-1) + uice(i+1,j-1) + - uice(i,j-1)) + + etaB(i+1,j) * ( uice(i+1,j) + vice(i+1,j) + - uice(i+1,j-1) ) + + etaB(i,j) * ( uice(i,j-1) + vice(i-1,j) ) ) + / Deltax ** 2 R2pp(i,j) = R2p(i,j) - vterm(i,j) * Sumu / 4.0 + R2pptmp c------------------------------------------------------------------------ c Calculate the x-comp of velocity: u_i,j c------------------------------------------------------------------------ case = maskB(i,j) + maskB(i,j+1) if ( case .eq. 0 ) then ! land uice(i,j) = 0.0 elseif ( case .eq. 2 ) then ! ice or water uice(i,j) = w * uice(i,j) + ( 1. - w ) + * ( A11(i,j) * R1pp(i,j) + A12(i,j) * R2pp(i,j) ) + / detA(i,j) elseif ( case .eq. 1 ) then ! ice or water adjacent to land deno = CdwC1(i,j) * costheta_w + + ( etaC(i,j) + etaC(i-1,j) + etaB(i,j+1) + etaB(i,j) ) + / deltax ** 2 RHStmp = ( etaC(i,j) * ( uice(i+1,j) - vice(i,j+1) + + vice(i,j) ) + + etaC(i-1,j) * ( uice(i-1,j) + vice(i-1,j+1) + - vice(i-1,j) ) + + etaB(i,j+1) * ( vice(i,j+1) - vice(i-1,j+1) ) + + etaB(i,j) * ( vice(i-1,j) - vice(i,j) ) ) + / deltax ** 2 if ( maskB(i,j) .eq. 0 ) then ! land just south RHStmp = RHStmp + + etaB(i,j+1) * uice(i,j+1) / Deltax ** 2 + + etaB(i,j) * uice(i,j+1) / 3.0 / Deltax ** 2 deno = deno + 2.0 * etaB(i,j) / Deltax ** 2 elseif ( maskB(i,j+1) .eq. 0 ) then ! land just north RHStmp = RHStmp + + etaB(i,j+1) * uice(i,j-1) / 3.0 / Deltax ** 2 + + etaB(i,j) * uice(i,j-1) / Deltax ** 2 deno = deno + 2.0 * etaB(i,j+1) / Deltax ** 2 endif uice(i,j) = w * uice(i,j) + ( 1. - w ) + * ( R1p(i,j) + uterm(i,j) * ( vice(i,j) + Sumv ) + / 4.0 + RHStmp ) / deno endif c------------------------------------------------------------------------ c Calculate the y-comp of velocity: v_i,j c------------------------------------------------------------------------ case = maskB(i,j) + maskB(i+1,j) if ( case .eq. 0 ) then ! land vice(i,j) = 0.0 elseif ( case .eq. 2 ) then ! ice or water vice(i,j) = w * vice(i,j) + ( 1. - w ) + * ( A21(i,j) * R1pp(i,j) + A22(i,j) * R2pp(i,j) ) + / detA(i,j) elseif ( case .eq. 1 ) then ! ice or water adjacent to land deno = CdwC2(i,j) * costheta_w + + ( etaB(i+1,j) + etaB(i,j) + etaC(i,j) + etaC(i,j-1) ) + / deltax ** 2 RHStmp = ( etaB(i+1,j) * ( uice(i+1,j) - uice(i+1,j-1) ) + + etaB(i,j) * (-uice(i,j) + uice(i,j-1) ) + + etaC(i,j) * ( vice(i,j+1) - uice(i+1,j) + + uice(i,j) ) + + etaC(i,j-1) * ( vice(i,j-1) + uice(i+1,j-1) + - uice(i,j-1) ) ) / Deltax ** 2 if ( maskB(i,j) .eq. 0 ) then ! land just east deno = deno + 2.0 * etaB(i,j) / Deltax ** 2 RHStmp = RHStmp + + etaB(i+1,j) * vice(i+1,j) / Deltax ** 2 + + etaB(i,j) * vice(i+1,j) / Deltax ** 2 / 3.0 elseif ( maskB(i+1,j) .eq. 0 ) then ! land just west deno = deno + 2.0 * etaB(i+1,j) / Deltax ** 2 RHStmp = RHStmp + + etaB(i+1,j) * vice(i-1,j) / Deltax ** 2 / 3.0 + + etaB(i,j) * vice(i-1,j) / Deltax ** 2 endif vice(i,j) = w * vice(i,j) + (1. - w) + * ( R2p(i,j) - vterm(i,j) * ( uice(i,j) + Sumu ) + / 4.0 + RHStmp ) / deno endif c------------------------------------------------------------------------ c Maximum error between succesive iterations c------------------------------------------------------------------------ if (abs(upast - uice(i,j)) .gt. error) then error = abs(upast - uice(i,j)) ierroru = i jerroru = j uflag = 1 vflag = 0 endif if (abs(vpast - vice(i,j)) .gt. error) then error = abs(vpast - vice(i,j)) ierrorv = i jerrorv = j uflag = 0 vflag = 1 endif enddo enddo k = k + 1 if (k .eq. nmax) then print *, 'maximum number of iteration exceeded in UVsolve' if ( uflag .eq. 1) + print *, 'max error on u = ', error, + 'at (i,j) =', ierror, jerror if ( vflag .eq. 1) + print *, 'max error on v = ', error, + 'at (i,j) =', ierror, jerror goto 10 endif if (error .gt. tol) goto 5 10 continue c print *, '# of iteration in UVsolve', k return end c**************************************************************************** c subroutine freedrift: calculates the freedrift solution c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c**************************************************************************** subroutine freedrift(nx,ny,ntrac) include 'comm_dyice.h' do i = 1, nx+1 do j = 1, ny+1 uice(i,j) = ( ufd(i,j) + ufd(i-1,j) ) / 2.0 + * min( 1, maskB(i,j) + maskB(i,j+1) ) vice(i,j) = ( vfd(i,j) + vfd(i,j-1) ) / 2.0 + * min( 1, maskB(i,j) + maskB(i+1,j) ) enddo enddo return end c**************************************************************************** c subroutine cavitating: c corrects the pressure and velocity fields calculated in 'UVsolve' c according to the following constraint c / c | = 0, if 0 < p < P_max, c div( u_i ) < < 0, if p = P_max, c | > 0, if p = 0. c \ c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c**************************************************************************** subroutine cavitating(nx,ny,ntrac) include 'comm_dyice.h' tol = 1e-05 ! tolerance on 'p' nmax = 2000 ! max number of iterations k = 0 ! counter c------------------------------------------------------------------------ c Maximum allowable internal ice pressure for a given 'h' and 'A' c------------------------------------------------------------------------ do i = 1, nx do j = 1, ny Pmax(i,j) = gradp * h(i,j) * exp(-C * ( 1 - A(i,j) ) ) enddo enddo 10 k = k+1 error = 0. do i = 1, nx do j = 1, ny div = ( ( uice(i+1,j) - uice(i,j) ) + + ( vice(i,j+1) - vice(i,j) ) ) / Deltax deno = ( 1. / CdwC1(i,j) + 1. / CdwC1(i+1,j) ) + / Deltax ** 2 + + ( 1. / CdwC2(i,j) + 1. / CdwC2(i,j+1) ) + / Deltay ** 2 pprime = -costheta_w * ( div - effshear(i,j) * tandelta ) + / deno c------------------------------------------------------------------------ c cap p to P_max (or p' to P_max - p) in convergence c cap p to 0 (or p' to 0 - p) in divergence c equation (13), Flato & Hibler (1992), equation ?? in c------------------------------------------------------------------------ pprime = min( pprime, Pmax(i,j) - p(i,j) ) + * max( int(-sign( 1., div ) ), 0 ) ! 0: div, 1: conv + + max( pprime, -p(i,j) ) + * max( int( sign( 1., div ) ), 0 ) ! 1: div, 0: conv c------------------------------------------------------------------------ c new ice velocity and internal ice pressure c boundary conditions on u_ice is set through the 'min' statement c------------------------------------------------------------------------ ufac = pprime / ( costheta_w * Deltax ) vfac = pprime / ( costheta_w * Deltay ) uice(i,j) = ( uice(i,j) - ufac / CdwC1(i,j) ) * + min( maskB(i,j) + maskB(i,j+1), 1 ) uice(i+1,j) = ( uice(i+1,j) + ufac / CdwC1(i+1,j) ) * + min( maskB(i+1,j) + maskB(i+1,j+1), 1 ) vice(i,j) = ( vice(i,j) - vfac / CdwC2(i,j) ) * + min( maskB(i,j) + maskB(i+1,j), 1 ) vice(i,j+1) = ( vice(i,j+1) + vfac / CdwC2(i,j+1) ) * + min( maskB(i,j+1) + maskB(i+1,j+1), 1 ) p(i,j) = p(i,j) + pprime c------------------------------------------------------------------------ c node where the maximum error occurs c------------------------------------------------------------------------ if (abs(pprime) .gt. error) then error = abs(pprime) ierror = i jerror = j endif enddo enddo if (k .eq. nmax) then print *, 'maximum number of iteration exceeded in cavitating' print *, 'max error = ', error, 'at (i,j) =', ierror, jerror goto 20 endif if (abs(error) .gt. tol) goto 10 20 continue c print *, '# of iteration in cavitating = ', k return end c************************************************************************ c Subroutine GranViscCoef: c computes the non-linear viscous coefficient at the nodes (etaB) c and the center (etaC) of the grid. c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c************************************************************************ subroutine GranViscCoef (nx,ny,ntrac,kk) include 'comm_dyice.h' c------------------------------------------------------------------------ c free slip boundary condition: c d(v_tangential)/d(normal) = 0 & v_normal = 0, at close boundary c------------------------------------------------------------------------ if ( BndyCond .eq. 2 ) then do i = 1, nx do j = 1, ny dudx = ( uice(i+1,j) - uice(i,j) ) / Deltax dvdy = ( vice(i,j+1) - vice(i,j) ) / Deltay dudyy = ( + ( uice(i,j+1) - uice(i,j) ) * maskB(i,j+1) + + ( uice(i+1,j+1) - uice(i+1,j) ) * maskB(i+1,j+1) + + ( uice(i,j) - uice(i,j-1) ) * maskB(i,j) + + ( uice(i+1,j) - uice(i+1,j-1) ) * maskB(i+1,j) ) + / (4.0 * Deltay) dvdxx = ( + ( vice(i+1,j) - vice(i,j) ) * maskB(i+1,j) + + ( vice(i+1,j+1) - vice(i,j+1) ) * maskB(i+1,j+1) + + ( vice(i,j) - vice(i-1,j) ) * maskB(i,j) + + ( vice(i,j+1) - vice(i-1,j+1) ) * maskB(i,j)) + / (4.0 * Deltay) deno = sqrt( ( dudx - dvdy ) ** 2 + + ( dudyy + dvdxx ) ** 2) deno = max( deno, 1e-15 ) etaC(i,j) = min ( p(i,j) * sinphi / deno, etamax ) enddo enddo do i = 2, nx do j = 2, ny etaB(i,j) = ( etaC(i,j) + etaC(i-1,j) + etaC(i,j-1) + + etaC(i-1,j-1) ) / 4.0 * maskB(i,j) enddo enddo c------------------------------------------------------------------------ c no-slip boundary condition: c v_normal = 0 and v_tangential = 0 c------------------------------------------------------------------------ elseif ( BndyCond .eq. 1) then do i = 1, nx+1 do j = 1, ny+1 dudx = 0.0 dvdy = 0.0 dudy(i,j) = 0.0 dvdx(i,j) = 0.0 if ( maskB(i,j) .eq. 1 ) then dudx = ( ( uice(i+1,j) + uice(i+1,j-1) ) - + ( uice(i-1,j) + uice(i-1,j-1) ) ) + / ( 4.0 * Deltax ) dvdx(i,j) = ( vice(i,j) - vice(i-1,j) ) / Deltax elseif ( maskB(i+1,j) - maskB(i-1,j) .eq. 1 ) then dudx = ( 4.0 * ( uice(i+1,j) + uice(i+1,j-1) ) + - ( uice(i+2,j) + uice(i+2,j-1) ) ) + / ( 4.0 * Deltax ) dvdx(i,j) = ( 9.0 * vice(i,j) - vice(i+1,j) ) + / ( 3.0 * Deltax ) elseif ( maskB(i+1,j) - maskB(i-1,j) .eq. -1 ) then dudx = (-4.0 * ( uice(i-1,j) + uice(i-1,j-1) ) + + ( uice(i-2,j) + uice(i-2,j-1) ) ) + / ( 4.0 * Deltax ) dvdx(i,j) = (-9.0 * vice(i-1,j) + vice(i-2,j) ) + / ( 3.0 * Deltax ) endif if (maskB(i,j) .eq. 1) then dvdy = ( ( vice(i-1,j+1) + vice(i,j+1) ) - + ( vice(i-1,j-1) + vice(i,j-1) ) ) + / ( 4.0 * Deltay ) dudy(i,j) = ( uice(i,j) - uice(i,j-1) ) / Deltay elseif ( maskB(i,j+1) - maskB(i,j-1) .eq. 1 ) then dvdy = ( 4.0 * ( vice(i-1,j+1) + vice(i,j+1) ) + - ( vice(i-1,j+2) + vice(i,j+2) ) ) + / ( 4.0 * Deltay ) dudy(i,j) = ( 9.0 * uice(i,j) - uice(i,j+1) ) + / ( 3.0 * Deltay ) elseif ( maskB(i,j+1) - maskB(i,j-1) .eq. -1 ) then dvdy = (-4.0 * ( vice(i-1,j-1) + vice(i,j-1) ) + + ( vice(i-1,j-2) + vice(i,j-2) ) ) + / ( 4.0 * Deltay ) dudy(i,j) = (-9.0 * uice(i,j-1) + uice(i,j-2) ) + / ( 3.0 * Deltay ) endif pice = ( p(i,j) + p(i-1,j) + p(i-1,j-1) + p(i,j-1) ) + / max( maskC(i,j) + maskC(i-1,j) + maskC(i-1,j-1) + + maskC(i,j-1), 1 ) deno = sqrt( ( dudx - dvdy ) ** 2 + + ( dudy(i,j) + dvdx(i,j) ) ** 2 ) deno = max( deno, 1e-15 ) etaB(i,j) = min ( pice * sinphi / deno, etamax ) enddo enddo do i = 1, nx do j = 1, ny dudx = ( uice(i+1,j) - uice(i,j) ) / Deltax dvdy = ( vice(i,j+1) - vice(i,j) ) / Deltax dudyc= ( dudy(i,j) + dudy(i+1,j) + dudy(i,j+1) + + dudy(i+1,j+1) ) / 4.0 dvdxc= ( dvdx(i,j) + dvdx(i+1,j) + dvdx(i,j+1) + + dvdx(i+1,j+1) ) / 4.0 deno = sqrt( ( dudx - dvdy ) ** 2 + ( dudyc + dvdxc ) ** 2 ) deno = max( deno, 1e-15 ) etaC(i,j) = min ( p(i,j) * sinphi / deno, etamax ) if ( kk .eq. 1 ) then effshear(i,j) = sqrt( ( dudx - dvdy ) ** 2 + + ( dudyc + dvdxc ) ** 2 ) / 2.0 endif enddo enddo endif return end c***************************************************************************** c Subroutne fluxx: c Computes the flux of a given tracer in the x-direction using the c upwind advection scheme c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c***************************************************************************** function fluxx (nx,ny,ntrac, + x, ! tracer being advected + i, ! ith coord of the grid cell + j, ! jth coord of the grid cell + k ! tracer number + ) include 'comm_dyice.h' dimension x(0:nx+1,0:ny+1,10) c------------------------------------------------------------------------ c x-coordinate of the grid cell upwind (ii) c------------------------------------------------------------------------ ii = i - max ( 0, int( sign( 1., uice(i,j) ) ) ) c------------------------------------------------------------------------ c check if the CFL criteria is met c------------------------------------------------------------------------ if ( uice(i,j) .gt. Deltax/Deltat ) + print *, 'WARNING: uice > Deltax/Deltat' c------------------------------------------------------------------------ c flux calculation c------------------------------------------------------------------------ fluxx = uice(i,j) * x(ii,j,k) return end c**************************************************************************** c Subroutine fluxy: c Compute the flux of a given tracer in the y-direction using the c upwind advection scheme c c Revision History c ---------------- c c Ver Date (dd-mm-yy) Author c c V01 14-05-97 L.-B. Tremblay c first version c c Address : Lamont-Doherty Earth Observatory of Columbia University c ------- Palisades, NY 10964-8000 USA c Email : tremblay@ldeo.columbia.edu c c**************************************************************************** function fluxy (nx,ny,ntrac, + x, ! tracer being advected + i, ! ith coord of the grid cell + j, ! jth coord of the grid cell + k ! tracer number + ) include 'comm_dyice.h' dimension x(0:nx+1,0:ny+1,10) c------------------------------------------------------------------------ c y-coordinate of the grid cell (jj) c------------------------------------------------------------------------ jj = j - max ( 0, int( sign( 1., vice(i,j) ) ) ) c------------------------------------------------------------------------ c check if the CFL criteria is met c------------------------------------------------------------------------ if ( vice(i,j) .gt. Deltax/Deltat ) + print *, 'WARNING: vice > Deltax/Deltat' c------------------------------------------------------------------------ c flux calculation c------------------------------------------------------------------------ fluxy = vice(i,j) * x(i,jj,k) return end c**************************************************************************** subroutine lsm_init(nx,ny,npt,iox,lsm) c**************************************************************************** dimension iox(npt),lsm(nx*ny) do i = 1, npt lsm(iox(i)) = i enddo return end