c$Source: /usr/our/senya/work/model/MC_PG/senq/RCS/dyn_new.f,v $ c$Author: senya $ c$Revision: 0.4 $ c$Date: 94/01/24 11:04:47 $ c$State: Exp $ c idig.f / "uphi" - model/ c-------------------------------------- function idig(xxxxx, n) integer xxxxx, n, res integer idig c res = xxxxx if (n .gt. 1) res = xxxxx / 10**(n-1) idig = mod (res, 10) return end function mlen (string) c--------------------------- character*(*) string integer mlen data LMAX /80/ k = 1 do while ( * k .le. LMAX .and. * int(string(k:k)) .ge. 33 .and. * int(string(k:k)) .le. 126) k = k + 1 enddo mlen = k-1 return end c ------------------------------------------------------------------ subroutine inout (ioin) c ------------------------------------------------------------------ c unix routine for opening the input parameter file and the output c logfil from the command line arguments. include 'comm_new.h' character*80 arg narg = iargc() i = 0 in = 0 do while (i .lt. narg) i = i + 1 call getarg (i, arg) if (arg .eq. '-h' .or. arg .eq. '-help') then goto 100 elseif (arg .eq. '-i') then i = i + 1 call getarg (i, finp) open (unit = ioin, file = finp) fout = finp(1:mlen(finp))//'.log' fcpu = finp(1:mlen(finp))//'.cpu' ftios = '.tios\0' in = 1 elseif (arg .eq. '-o') then i = i + 1 call getarg (i, fout) elseif (arg .eq. '-t') then i = i + 1 call getarg (i, ftios) elseif (arg .eq. '-d') then i = i + 1 call getarg (i, fbi) n_in = mlen(fbi) call dump_rstrt stop endif enddo if (in .eq. 0) goto 100 open (unit = iout, file = fout) goto 200 100 call getarg (0, arg) write (6, *) * 'usage: '//arg(1:mlen(arg))//' [-i file][-t ftios][-d file]' write (6,*) 'where: -i file - for model control ' write (6,*) ' -t file - for tios control (deflt:<.tios>)' write (6,*) ' -d file - make a dump of data/restart ' stop 200 return end c ------------------------------------------------------------------ subroutine ddiv (npt,nzi,uc,vc,emx,emy,emxy,w,fhd,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,tp,mbc,lpbcwk,lpbcek,h,bdiv) csenq ------------------------------------------------------------------ c compute the divergence (fhd) for all layers and put the Sum in w(1,nz). implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz), * snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz), * lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz) dimension uc(npt,nz),vc(npt,nz),emx(npt),emy(npt),emxy(npt),w(npt,nz), * fhd(npt,nz), tp(npt,4),h(npt,nz), bdiv(npt), nzi(npt) c c set boundary condition flag based on whether interior corners are c treated as boundaries. see bcset and dfdx. c nbu = 0 nbv = 0 if(mbc.eq.1 .or. mbc.eq.4) nbu = 1 if(mbc.eq.1 .or. mbc.eq.3) nbv = 1 c........compute d(hv)/dy & d(hu)/dx.. nxk = nbxk(1) nyk = nbyk(1) nck = ncsk(1) npbk = npbck(1) call dfdx1(uc,tp(1,3),npt,nbu,nxk,nyk,nck,lxxk,lyxk, * snxk,npbk,lpbcwk,lpbcek) call dfdy1(vc,tp(1,4),npt,nbv,nyk,nxk,nck, * lyyk,lxyk,snyk,isyk) if (mgrid .ne. 2) then do i = 1, npt fhd(i,1) = emx(i)*tp(i,3) + emy(i)*tp(i,4) enddo else do i = 1, npt fhd(i,1) = emx(i)*tp(i,3) + emy(i)*tp(i,4) + emxy(i)*vc(i,1) enddo endif do k = 2, nz npk = nptk(k) c........mud points have zero transport: call zero_em (npt, tp) call zero_em (npt, tp(1,2)) do j = 1, npk i = isk(j,k) tp(i,1) = uc(i,k) tp(i,2) = vc(i,k) enddo call dfdx1(tp,tp(1,3),npt,nbu,nxk,nyk,nck,lxxk,lyxk, * snxk,npbk,lpbcwk,lpbcek) call dfdy1(tp(1,2),tp(1,4),npt,nbv,nyk,nxk,nck, * lyyk,lxyk,snyk,isyk) c........now multiply by the appropriate scale factors to find divergence. c........div(u) = (1/mx)*(du/dx) + (1/my)*(dv/dy) + myx*u + mxy*v c........we also accumulate the sum of layer divergences in w(nz) if (mgrid .ne. 2) then do i = 1, npt fhd(i,k) = emx(i)*tp(i,3) + emy(i)*tp(i,4) enddo else do i = 1, npt fhd(i,k) = emx(i)*tp(i,3) + emy(i)*tp(i,4) + * emxy(i)*tp(i,2) enddo endif enddo do i = 1, npt mz = nzi(i) fhd(i,1) = fhd(i,1) - h(i,1)*bdiv(i) w(i,1) = fhd(i,1) do k = 2, mz fhd(i,k) = fhd(i,k) - h(i,k)*bdiv(i) c w(i,k) = w(i,k-1) + fhd(i,k) enddo enddo c do k = 1, nz c do i = 1, npt c w(i,k) = 0. c enddo c enddo return end c ---------------------------------------------------------- subroutine wtop (npt, w, fhd, fh, h, ncyc, istep, dnt) c ---------------------------------------------------------- include 'comm_new.h' c.....subroutine to set W at first interface. dimension fhd(1), w(1), fh(1), h(1) if (imix .eq. 4) then if (mix_wtop .eq. 1) then do i = 1, npt w(i) = 0. enddo elseif (mix_wtop .eq. 2) then binv1 = real(ncyc-istep)/dnt do i = 1, npt dmin = fh(i) + binv1 * (h(i) - hmax_mix) dmax = fh(i) + binv1 * (h(i) - hmin_mix) divw = fhd(i) divw = amin1(divw, dmax) divw = amax1(divw, dmin) w(i) = fhd(i) - divw enddo endif else if (iv_top .eq. 1) then do i = 1, npt w(i) = fhd(i) enddo elseif (iv_top.eq.3) then do i = 1, npt w(i) = 0. enddo endif endif return end c ------------------------------------------------------------------ subroutine dwcal (npt, nz, nsig, nzi, w, fhd, sigma) c ------------------------------------------------------------------ c.....find w(k+1/2) at each intermideate sigma-layer interface c.....uses sum of div() from w(nz) implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension w(npt,nz),fhd(npt,nz),sigma(nz),nzi(npt) do i = 1, npt win = - fhd(i,1) + w(i,1) do k = 2, nzi(i) - 1 w(i,k) = w(i,k-1) + fhd(i,k) - 1.5*(sigma(k)+sigma(k+1))*win enddo c to check that w(k=last) = 0. c k = nzi(i) c w(i,k) = w(i,k-1) + fhd(i,k) - 1.5*sigma(k)*win enddo return end c ------------------------------------------------------------------ subroutine vertu (npt,nz,nsig,nzi,nzi_b,bint,taux,tauy,u,v,w,h, * fu,fv,fh,vertx,verty,zfu,zfv) c----------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension bint(1), taux(npt), tauy(npt), u(npt,nz),v(npt,nz), * w(npt,nz), h(npt,nz), fu(npt,nz), fv(npt,nz), nzi(npt) * , zfu(npt), zfv(npt), nzi_b(npt), fh(npt,nz) #ifdef dump_all * , vertx(npt,nz), verty(npt,nz) #endif c.....surface layer: #ifdef dump_all do i = 1, npt do k = 1, nz vertx(i,k) = 0. verty(i,k) = 0. enddo enddo #endif do i = 1, npt fu(i,1) = fu(i,1) + taux(i) fv(i,1) = fv(i,1) + tauy(i) enddo do k = 1, nsig k1 = k + 1 do i = 1, npt w_add = w(i,k) fh(i,k) = fh(i,k) + w_add fh(i,k1) = fh(i,k1) - w_add enddo enddo c.....going by layer interfaces: do i = 1, npt mz = nzi(i) dz = h(i,1) do k = 1, mz - 1 b2 = bint(k) k1 = k + 1 dz = 2.*h(i,k) - dz h_ave = b2/dz u_ave = 0.5*(u(i,k) + u(i,k1)) v_ave = 0.5*(v(i,k) + v(i,k1)) w_add = w(i,k) u_add = h_ave*(u(i,k1) - u(i,k)) + w_add*u_ave v_add = h_ave*(v(i,k1) - v(i,k)) + w_add*v_ave fu(i,k) = fu(i,k) + u_add fu(i,k1) = fu(i,k1) - u_add fv(i,k) = fv(i,k) + v_add fv(i,k1) = fv(i,k1) - v_add if (k.eq.nzi_b(i)) then zfu(i) = zfu(i) + u_add zfv(i) = zfv(i) + v_add endif #ifdef dump_all vertx(i,k) = vertx(i,k) + u_add vertx(i,k1) = vertx(i,k1) + -u_add verty(i,k) = verty(i,k) + v_add verty(i,k1) = verty(i,k1) - v_add #endif enddo enddo return end subroutine vertt (npt,nz,nzi,cint,q,qr,w,h,t,ft) c---------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension cint(nz),q(npt),qr(npt),w(npt,nz),h(npt,nz), * t(npt,nz),ft(npt,nz),tp(1) dimension nzi(npt) parameter (H_ATT1 = -1./17.) do i = 1, npt ft(i,1) = ft(i,1) + q(i) + qr(i) dz = h(i,1) do k = 1, nzi(i) - 1 dz = 2.*h(i,k) - dz k1 = k + 1 c2 = cint(k) tk = t(i,k) tk1 = t(i,k1) h_ave = c2/dz t_ave = 0.5*(tk + tk1) t_add = h_ave*(tk1 - tk) + w(i,k)*t_ave ft(i,k) = ft(i,k) + t_add ft(i,k1) = ft(i,k1) - t_add enddo enddo if (isolrp .eq. 1) then c.....add penetrating solar radiation: do i = 1, npt z = 0. do k = 1, nzi(i) - 1 k1 = k + 1 z = z + h(i,k) t_add = qr(i)*exp(H_ATT1*z) ft(i,k) = ft(i,k) - t_add ft(i,k1) = ft(i,k1) + t_add enddo enddo endif return end subroutine vertts (npt,nz,nzi,cint,q,qr,ep,w,h,t,ft,s,fs) c------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension cint(nz), q(npt), qr(npt), ep(npt), w(npt,nz), h(npt,nz), * t(npt,nz), ft(npt,nz), s(npt,nz), fs(npt,nz), tp(1) dimension nzi(npt) parameter (H_ATT1 = -1./17.) do i = 1, npt ft(i,1) = ft(i,1) + q(i) + qr(i) fs(i,1) = fs(i,1) + ep(i) dz = h(i,1) do k = 1, nzi(i) - 1 dz = 2.*h(i,k) - dz k1 = k + 1 c2 = cint(k) h_ave = c2/dz tk = t(i,k) tk1 = t(i,k1) t_ave = 0.5*(tk + tk1) sk = s(i,k) sk1 = s(i,k1) s_ave = 0.5*(sk + sk1) wik = w(i,k) t_add = h_ave*(tk1 - tk) + wik*t_ave s_add = h_ave*(sk1 - sk) + wik*s_ave ft(i,k) = ft(i,k) + t_add ft(i,k1) = ft(i,k1) - t_add fs(i,k) = fs(i,k) + s_add fs(i,k1) = fs(i,k1) - s_add enddo enddo if (isolrp .eq. 1) then c.....add the penetrating solar radiation: do i = 1, npt z = 0. do k = 1, nzi(i) - 1 k1 = k + 1 z = z + h(i,k) t_add = qr(i)*exp(H_ATT1*z) ft(i,k) = ft(i,k) - t_add ft(i,k1) = ft(i,k1) + t_add enddo enddo endif return end subroutine verttr (npt,nz,nzi,cint,w,h,tr,ftr) c------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension cint(nz), w(npt,nz), h(npt,nz), tr(npt,nz,1), ftr(npt,nz,1) dimension nzi(npt) parameter (H_ATT1 = -1./17.) do i = 1, npt dz = h(i,1) do k = 1, nzi(i) - 1 dz = 2.*h(i,k) - dz k1 = k + 1 c2 = cint(k)/2. h_ave = c2/dz wik = w(i,k)/2. do n = 1, ntrac trk = tr(i,k,n) trk1 = tr(i,k1,n) tr_ave = trk + trk1 tr_add = h_ave*(trk1 - trk) + wik*tr_ave ftr(i,k,n) = ftr(i,k,n) + tr_add ftr(i,k1,n) = ftr(i,k1,n) - tr_add enddo enddo enddo return end subroutine enso2date (enso, id, im, iy) c---------------------------------------------------- #define LEAP_YEAR(y) (mod(y,4) .eq. 0) integer*2 norm(12) data norm /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ iy = 1960 + jint(enso/12.) res = enso - (iy - 1960.)*12. im = jint(res) + 1 if (im .eq. 2 .and. LEAP_YEAR(iy) ) then id = 1 + int(29 * (res - int(res))) else id = 1 + int(norm(im) * (res - int(res))) endif return end function date2enso (id, im, iy) c-------------------------------------- #define LEAP_YEAR(y) (mod(y,4) .eq. 0) integer*2 norm(12) data norm /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ if (im .eq. 2 .and. LEAP_YEAR(iy) ) then date2enso = (iy - 1960.) * 12. + 1. + real((id-1))/29. else date2enso = (iy - 1960.) * 12. + real(im-1) + real((id-1))/norm(im) endif return end subroutine DayOfYear(enso, idoy, idiy) c------------------------------------------------ #define LEAP_YEAR(y) (mod(y,4) .eq. 0) integer*2 norm(12) data norm /0, 31, 59, 90,120,151,181,212,243,273,304,334/ call enso2date(enso, id, im, iy) if (LEAP_YEAR(iy)) then idiy = 366 else idiy = 365 endif if (im .gt. 2 .and. LEAP_YEAR(iy) ) then idoy = norm(im) + id + 1 else idoy = norm(im) + id endif return end subroutine enso2res (renso, id, im, iy) c---------------------------------------------------- #define LEAP_YEAR(y) (mod(y,4) .eq. 0) integer*2 norm(12) data norm /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ iy = jint(renso/12.) res = renso - real(12*iy) im = jint(res) res = abs(res - real(im)) id = jint(res*norm(im+1)) return end c ------------------------------------------------------------------ subroutine knergy(npt,nz,nptk,isk,area,basin,h,u,v,en) c ------------------------------------------------------------------ c subroutine to compute the kinetic energy c c the kinetic energy is given by c k.e. = sum < 1/2 h(k)*u(k)^2 >. c c en(1,k) = (output) kinetic energy for layer k. c en(2,k) = not used c note: area(i) = .5*dx*dy c note: basin = sum(area) (i.e., half the surface area) c implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension u(npt,1),v(npt,1),h(npt,1),area(1),en(npten,1), * nptk(1),isk(npt,1) c........Compute the Kinetic Energy, each layer and total eksum = 0.0 do k = 1, nz ek = 0.0 do j = 1, nptk(k) i = isk(j,k) uu = u(i,k) vv = v(i,k) ek = ek + area(i) * h(i,k) * (uu*uu + vv*vv) enddo en(1,k)= 0.5 * ek / basin eksum = eksum + en(1,k) enddo en(1,nz+1) = eksum return end c ------------------------------------------------------------------ subroutine pnergy (ifrst,npt,nptk,nz,isk,h,area,t,dens,w,basin,en,tp) c ------------------------------------------------------------------ c subroutine to compute the potential energy, the heat content, c and the temperature variance and total mass. c en(k,3) = (output) Potential Energy for layer k. c en(k,4) = (output) Heat Content for layer k. c en(k,5) = (output) Mass Content for layer k. c en(k,6) = (output) Volume for layer k. c p.e. = sum - 1/2 g , c b(k) = alpha*g*(t(k)-t(b)) c Heat Content = sum c Mass Content = sum c Volume = sum c include 'comm_para.h' include 'comm_new.h' logical ifrst dimension h(npt,1),area(1),t(npt,1),dens(npt,1),w(npt,1), * en(npten,1), tp(1), tmp(npt,1), nptk(1), isk(npt,1) pointer (ptmp, tmp) if (use_salt) then ptmp = loc(dens) base = SITUD_BOT shift = 1000. else ptmp = loc(t) base = TEMP_BOT shift = 1. endif do i = 1, npt tp(i) = h(i,1) enddo epsum = 0. hcsum = 0. wcsum = 0. vlsum = 0. c.......................integrate over each layer. do k = 1, nz epk = 0. hck = 0. wck = 0. vlk = 0. do j = 1, nptk(k) i = isk(j,k) hk = h(i,k) ahk = area(i) * hk ahbk = ahk * (tmp(i,k) - base) epk = epk + ahbk * (hk - w(i,nz) + tp(i)) hck = hck + ahbk wck = wck + ahk * (shift + tmp(i,k)) vlk = vlk + ahk enddo en(3,k) = epk en(4,k) = hck en(5,k) = wck en(6,k) = vlk epsum = epsum + epk hcsum = hcsum + hck wcsum = wcsum + wck vlsum = vlsum + vlk c find 2*(sum h(j)+h(k)/2) for the next layer. if (k .lt. nz) then do j = 1, nptk(k) i = isk(j,k) tp(i) = tp(i) + h(i,k) + h(i,k+1) enddo endif enddo if ( ifrst ) then epf1 = 1./epsum hcf1 = 1./hcsum wcf1 = 1./wcsum vlf1 = 1./vlsum endif do k = 1, nz en(3,k) = epf1 * en(3,k) en(4,k) = hcf1 * en(4,k) en(5,k) = wcf1 * en(5,k) en(6,k) = vlf1 * en(6,k) enddo en(3,nz+1) = epf1 * epsum en(4,nz+1) = hcf1 * hcsum en(5,nz+1) = wcf1 * wcsum en(6,nz+1) = vlf1 * vlsum return end c subroutine vel_updat(npt,nz,nzi,binv,abinv,uc,vc,fu,fv) c ---------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension uc(npt,nz),vc(npt,nz),fu(npt,nz),fv(npt,nz),nzi(npt) do i = 1, npt do k = 1, nzi(i) uc(i,k) = uc(i,k) + binv*fu(i,k) vc(i,k) = vc(i,k) + binv*fv(i,k) fu(i,k) = abinv*fu(i,k) fv(i,k) = abinv*fv(i,k) enddo enddo return end c ------------------------------------------------------------------ subroutine dhoriz(npt,u,v,uc,vc,f,fu,fv,fhd,emx,emy,emxy,tp,mbc,zfu,zfv, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,lpbcwk,lpbcek,nzi * ,corx,cory,xnl,ynl,fh) c ------------------------------------------------------------------ c subroutine that calculates the horizontal terms in the momentum c equation. e.g. the coriolis terms and the horizontal advection c terms. c c npt = (input) # of grid points per model layer (nxyc or nxy). c u,v = (input) zonal & merid. velocity for time step n. c uc,vc = (input) mass transport. c fu,fv = (input/output) update transport arrays c emx,emy= (input) factor for x,y-differencing, d(psix)/dx*1/delx. c lxx,...= (input) nbx+ncs indices of the ocean x,y-segment end points c snx,sny= (input) nbx+ncs signs (+1 or -1) c isy = (input) indices to convert from an x-sort to a y-sort. c tp = (input) temporary space. c mbc = (input) type of boundary condition on u and v c implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch dimension u(npt,nz),v(npt,nz),uc(npt,nz),vc(npt,nz),f(npt),fu(npt,nz),fv(npt,nz), * emx(npt),emy(npt),emxy(npt), tp(npt,1), zfu(npt), zfv(npt), fhd(npt,nz) * , corx(npt,nz), cory(npt,nz) * , xnl(npt,nz), ynl(npt,nz), fh(npt,nz) dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz), * snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz), * lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz),nzi(npt) c c set boundary condition flag for differencing at or near the c boundary. c c for cases where the derivative is in the direction to the flow c for the cases where u slips along zonal boundaries (mbc=2,3), c interior corners are not used. c nbu = 0 c c for no-slip everywhere (mbc=1) or at zonal boundaries (mbc=0), c interior corners are used for a one-sided difference. c if(mbc.eq.1 .or. mbc.eq.4) nbu = 1 c c for the cases where v slips along meridional boundaries (mbc=2,4), c interior corners are not used. c nbv = 0 c c for no-slip everywhere (mbc=1) or at meridional boundaries c (mbc=3), interior corners are used. c c.....update fh in sigma layers: do k = 1, nsig do i = 1, npt fh(i,k) = fh(i,k) - fhd(i,k) enddo enddo if(mbc.eq.1 .or. mbc.eq.3) nbv = 1 do k = 1, nz npk = nptk(k) nxk = nbxk(k) nyk = nbyk(k) nck = ncsk(k) npbk = npbck(k) c.....add coriolis terms: do j = 1, npk i = isk(j,k) tmp = f(i)*vc(i,k) fu(i,k) = fu(i,k) + tmp corx(i,k) = tmp tmp = -f(i)*uc(i,k) fv(i,k) = fv(i,k) + tmp cory(i,k) = tmp enddo c.....add d(hu^2)/dx do j = 1,npk i = isk(j,k) tp(i,2) = uc(i,k)*u(i,k) ynl(i,k) = 0. enddo call dfdxk(tp(1,2),tp,npt,npk,nbu,nxk,nyk,nck,lxxk(1,k),lyxk(1,k), * snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k)) do j = 1, npk i = isk(j,k) tmp = emx(i)*tp(i,1) fu(i,k) = fu(i,k) - tmp xnl(i,k) = -tmp enddo if (mgrid .eq. 2) then do j = 1, npk i = isk(j,k) tmp = emxy(i)*tp(i,2) fv(i,k) = fv(i,k) + tmp ynl(i,k) = +tmp enddo endif c.....find dv**2h/dy do j = 1, npk i = isk(j,k) tp(i,2) = vc(i,k)*v(i,k) enddo call dfdyk(tp(1,2),tp,npt,npk,nbv,nyk,nxk,nck,lyyk(1,k),lxyk(1,k), * snyk(1,k),isyk(1,k)) do j = 1, npk i = isk(j,k) tmp = emy(i)*tp(i,1) fv(i,k) = fv(i,k) - tmp ynl(i,k) = ynl(i,k)-tmp enddo if (mgrid .eq. 2) then do j = 1, npk i = isk(j,k) tmp = emxy(i)*tp(i,2) fv(i,k) = fv(i,k) - tmp ynl(i,k) = ynl(i,k)-tmp enddo endif c.....find d(u*h*v)/dx do j = 1, npk i = isk(j,k) tp(i,3) = uc(i,k)*v(i,k) enddo call dfdxk(tp(1,3),tp,npt,npk,1,nxk,nyk,nck,lxxk(1,k),lyxk(1,k), * snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k)) do j = 1, npk i = isk(j,k) tmp = emx(i)*tp(i,1) fv(i,k) = fv(i,k) - tmp ynl(i,k) = ynl(i,k)-tmp enddo c.....find d(u*h*v)/dy call dfdyk(tp(1,3),tp,npt,npk,1,nyk,nxk,nck,lyyk(1,k),lxyk(1,k), * snyk(1,k),isyk(1,k)) do j = 1, npk i = isk(j,k) tmp = emy(i)*tp(i,1) fu(i,k) = fu(i,k) - tmp xnl(i,k) = xnl(i,k)-tmp enddo if (mgrid .eq. 2) then do j = 1, npk i = isk(j,k) tmp = 2.*emxy(i)*tp(i,3) fu(i,k) = fu(i,k) - tmp xnl(i,k) = xnl(i,k)-tmp enddo endif enddo do i = 1, npt mz = nzi(i) do k = 1, mz c uncomment next 2 lines to take out the nonlinear terms in the barotropic part c zfu(i) = zfu(i) + corx(i,k) + xnl(i,k) c zfv(i) = zfv(i) + cory(i,k) + ynl(i,k) zfu(i) = zfu(i) + corx(i,k) zfv(i) = zfv(i) + cory(i,k) enddo enddo return c end of dhoriz. end c ------------------------------------------------------------------ subroutine thoriz(npt,uc,vc,t,ft,fhd,emx,emy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,lok,tp,mbc,lpbcwk,lpbcek) c ------------------------------------------------------------------ c note: mtc is not used in this version of the code. rather c n.grad(t)=0 is set at all closed boundaries. at open boundaries c determined in the input grid file (see rdgrid), the temperature c derivative is not zeroed and t is set in tbcset. c c subroutine that calculates the horizontal terms in the c temperature equation. c implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) dimension uc(npt,nz),vc(npt,nz),t(npt,nz),ft(npt,nz),emx(npt),emy(npt), * tp(npt,3),fhd(npt,nz) dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz), * snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz), * lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz),lok(4*MAXSID,nz) c c in order to impose the zero heat flux condition at the boundaries c the divergence operator is broken up. c c set boundary condition flag for differencing the mass transport c at or near the boundary in the direction of the flow. c c for the cases where u slips along zonal boundaries (mbc=2,3), c interior corners are not used. c nbu = 0 c c for no slip everywhere (mbc=1) or at zonal boundaries (mbc=4), c interior corners are used for a one-sided difference. c if(mbc.eq.1 .or. mbc.eq.4) nbu = 1 c c for the cases where v slips along meridional boundaries c (mbc=2,4), interior corners are not used. c nbv = 0 c c for no slip everywhere (mbc=1) or at meridiional boundaries c (mbc=3), interior corners are used. c if(mbc.eq.1 .or. mbc.eq.3) nbv = 1 c do k = 1, nz npk = nptk(k) nxk = nbxk(k) nyk = nbyk(k) nck = ncsk(k) npbk = npbck(k) nok = nlok(k) c.....compute dt/dx (dt/dx=0 at x-boundaries) call dfdxk(t(1,k),tp,npt,npk,0,nxk,nyk,nck,lxxk(1,k),lyxk(1,k), * snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k)) call zerodt(tp,nok,lok(1,k),nxk,lxxk(1,k),tp(1,3)) c.....compute dt/dy (dt/dy=0 at y-boundaries) call dfdyk(t(1,k),tp(1,2),npt,npk,0,nyk,nxk,nck,lyyk(1,k),lxyk(1,k), * snyk(1,k),isyk(1,k)) call zerodt(tp(1,2),nok,lok(1,k),nyk,lyxk(1,k),tp(1,3)) c.....update the heat content array. do j = 1, npk i = isk(j,k) ft(i,k) = ft(i,k)-(emx(i)*uc(i,k)*tp(i,1) + emy(i)*vc(i,k)*tp(i,2)) ft(i,k) = ft(i,k) - t(i,k)*fhd(i,k) enddo enddo return c end of thoriz. end c subroutine zero_em (n, a) c-------------------------------- dimension a(1) do i = 1, n a(i) = 0.0 enddo return end subroutine copya2b (n, a, b) c-------------------------------- dimension a(1), b(1) do i = 1, n b(i) = a(i) enddo return end c ------------------------------------ subroutine decap(npt, nz, nzi, u,v,uc,vc,h) c ------------------------------------ c subroutine to convert horizontal transport to horizontal velocity. implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension u(npt,nz),v(npt,nz),uc(npt,nz),vc(npt,nz),h(npt,nz),nzi(npt) c do i = 1, npt do k = 1, nzi(i) hinv = 1./h(i,k) u(i,k) = uc(i,k)*hinv v(i,k) = vc(i,k)*hinv enddo enddo return end c ------------------------------------------------------------------ subroutine capfrm(npt,nz,nzi, u,v,uc,vc,h) c ------------------------------------------------------------------ c subroutine to convert velocities to transport. c implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension u(npt,nz),v(npt,nz),uc(npt,nz),vc(npt,nz),h(npt,nz),nzi(npt) c do i = 1, npt do k = 1, nzi(i) hk = h(i,k) uc(i,k) = u(i,k)*hk vc(i,k) = v(i,k)*hk enddo enddo return c end of capfrm. end c ------------------------------------- subroutine tdecap(npt, nzi, t, h) c ------------------------------------- c.....convert heat content to temperature. implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension t(npt,1), h(npt,1), nzi(npt) do i = 1, npt do k = 1, nzi(i) t(i,k) = t(i,k)/h(i,k) enddo enddo return end subroutine fixed_dep (npt,nzi_b,h,fu,fv,tfu,tfv,rhsx,rhsy,crhsx,crhsy) c----------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension h(npt,1),fu(npt,1),fv(npt,1),tfu(1),tfv(1),nzi_b(npt) #ifdef dump_all dimension rhsx(npt,1), rhsy(npt,1) dimension crhsx(npt,1), crhsy(npt,1) #endif do i = 1, npt do k = 1, nzi_b(i) hk = h(i,k) #ifdef dump_all rhsx(i,k) = fu(i,k) rhsy(i,k) = fv(i,k) #endif fu(i,k) = fu(i,k) - hk* tfu(i) fv(i,k) = fv(i,k) - hk* tfv(i) #ifdef dump_all crhsx(i,k) = fu(i,k) crhsy(i,k) = fv(i,k) #endif enddo enddo return end c ------------------------------------------------------------------ subroutine capt(npt,nz,nzi,t,h) c ------------------------------------------------------------------ c form the heat content, ht, from the temperature. c dimension t(npt,nz),h(npt,nz),nzi(npt) c do i = 1, npt do k = 1, nzi(i) t(i,k) = t(i,k)*h(i,k) enddo enddo return c end of capt. end c ------------------------------------------------------------------ subroutine tupdat(npt,nz,nzi,binv,abinv,t,ft) c ------------------------------------------------------------------ c update temperature fields as was done for u, v, h in updat1. c implicit real(a-h,o-z),integer(i-n) dimension t(npt,nz),ft(npt,nz),nzi(npt) c do i = 1, npt do k = 1, nzi(i) t(i,k) = t(i,k) + binv*ft(i,k) ft(i,k) = abinv*ft(i,k) enddo enddo return c end of tupdat. end c ------------------------------------------------------------------ subroutine zerodt(dt,nlo,lo,nb,lb,tp) c ------------------------------------------------------------------ c set flux of t zero at the boundaries. c dt/dx = 0 at the x-sidewall boundaries. c dt/dy = 0 at the y-sidewall boundaries. c dimension dt(1),lo(1),tp(1),lb(1) c if (nlo.gt.0) then c first save the open boundary dt values. do i=1,nlo tp(i) = dt(lo(i)) enddo endif c c now zero dt at all x or y-boundaries. do i=1,nb dt(lb(i)) = 0. enddo c if (nlo.gt.0) then c replace dt values at open boundaries. do i=1,nlo dt(lo(i)) = tp(i) enddo endif c return c end of zerodt. end c ------------------------------------------------------------------ subroutine bcset(mbc,lxxk,lyyk,npt,u,v,nzi,nzi_b) c ------------------------------------------------------------------ c impose the u, v boundary conditions according to mbc. c c mbc = (input) type of boundary condition: c = 1; u(xb)=v(yb)=u(yb) = v(xb) = 0; no slip everywhere. c = 2; u(xb)=v(yb) = 0; no normal flow. c = 3; u(xb)=v(yb)=du(yb)/dy= v(xb) = 0; no slip at eastern c and western side walls; free slip along northern and c southern boundaries/steps, v=du/dy=0. c = 4; u(xb)=v(yb)=u(yb) = dv(xb)/dx= 0; no slip at northern c and southern; free slip along eastern and western c boundaries/steps, u=dv/dy=0. c c lxx = (input) nbx x-boundary plus ncs corner indices for a c regular or compressed x-sort. c lyy = (input) nby y-boundary plus ncs corner indices for a c regular or compressed x-sort. c npt = (input) number of field points/layer. c u,v = (input) fields. c = (output) fields with boundary conditions imposed. c implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) dimension lxxk(MXBDY,nz), lyyk(MXBDY,nz) dimension u(npt,nz),v(npt,nz),nzi(npt),nzi_b(npt) c c normal components are always zero at boundaries, u(xb)=v(yb)=0. c similarly, so is the along boundary derivative c du(xb)/dy = dv(yb)/dx = 0. c c the ncs corner points are part of the u/v-boundaries depending c on mbc. c mbc = 1; yes for u and v; du(yb)/dx = dv(xb)/dy = 0. c = 2; no for u and v. c = 3; no for u, yes for v; dv(xb)/dy = 0. c = 4; yes for u, no for v; du(yb)/dx = 0. c do k = 1, nz do i = 1, nbxk(k) u(lxxk(i,k),k) = 0. enddo do i = 1, nbyk(k) v(lyyk(i,k),k) = 0. enddo enddo do i = 1, npt do k = nzi_b(i)+1,nzi(i) u(i,k) = 0. v(i,k) = 0. enddo enddo if(mbc.eq.1 .or. mbc.eq.4) then do 30 k=1,nz do 30 i=1,nbyk(k)+ncsk(k) 30 u(lyyk(i,k),k) = 0. endif c if(mbc.eq.1 .or. mbc.eq.3) then do 40 k=1,nz do 40 i=1,nbxk(k)+ncsk(k) 40 v(lxxk(i,k),k) = 0. endif return c end of bcset. end c -------------------------------------------------------- logical function non_stable (iout, npt, nxp, nz, iox, t, u, v) c -------------------------------------------------------- c check to see that t or velocities are not bizarre c h = (input) layer thickness field. implicit real(a-h,o-z),integer(i-n) dimension t(npt,1), iox(1), u(npt,1), v(npt,1) c non_stable = .false. icheck = 0 do i = 1, npt if (t(i,2) .lt. -10..or.t(i,2).gt.50.or. * u(i,1)**2+v(i,1)**2.gt.400.) icheck = icheck + 1 enddo non_stable = (icheck .ne. 0) if (non_stable) then write (iout, *) 'Number of illegal points =', icheck do i = 1, npt if (t(i,2) .lt. -10..or.t(i,2).gt.50) then jj = 1 + (iox(i)-1)/nxp ii = iox(i) - (jj-1)*nxp write (iout, 11) i, ii, jj, t(i,2) endif if (u(i,1)**2 + v(i,1)**2 .gt. 400.) then jj = 1 + (iox(i)-1)/nxp ii = iox(i) - (jj-1)*nxp write (iout, 12) i, ii, jj endif enddo endif 11 format ('t(k=',i4,',2)[i=',i3,',j=',i3,'] =', g10.3) 12 format ('i,ii,jj=',3i8,' velocity is greater than 20 m/s') return end subroutine h_updat(npt,nsig,binv,abinv,h,fh) c------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension h(npt,1),fh(npt,1) include 'comm_new.h' do k = 1, nsig do i = 1, npt h(i,k) = h(i,k) + binv*fh(i,k) fh(i,k) = abinv*fh(i,k) enddo enddo return end c ------------------------------------------------------------------ subroutine btpgf(npt,nzi_b,h,temp,dens,fu,fv,emx,emy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,lok,tp,tq,tr,lpbcwk,lpbcek,zfu,zfv, * pgfx,pgfy) c ------------------------------------------------------------------ c subroutine that calculates the pressure gradient terms in the c momentum equation c implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' include 'comm_new.h' dimension h(npt,nz),temp(npt,nz), fu(npt,nz),fv(npt,nz), * emx(1),emy(1),tp(npt,4),tq(npt,4),tr(npt,1), * dens(npt,nz),tmp(npt,nz),zfu(npt),zfv(npt), * pgfx(npt,nz),pgfy(npt,nz) dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz), * snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),nzi_b(npt), * lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz),lok(4*MAXSID,nz) pointer (p_tmp, tmp) common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) c-new temp array should contain temp, not heat content!!! c-new dens array should contain dens, not mass!!! c if (use_salt) then c .......................................case of density from EOS: c b = (grav/(1000+sigma0)) * (sigma0 - sigma(k)) c coef = -GRAVTY / (1000. + SITUD_BOT) bottom = SITUD_BOT p_tmp = loc(dens) else c .....................case of linear (in T) density and buoyancy: c b = alph * grav * (t(k) - t_bot) coef = TALPHA * GRAVTY bottom = TEMP_BOT p_tmp = loc(temp) endif do i = 1, npt dh = h(i,1)/2. tp(i,1) = dh*tmp(i,1) tp(i,2) = dh tp(i,4) = dh tq(i,3) = coef * emx(i) tq(i,4) = coef * emy(i) enddo do k = 1, nz npk = nptk(k) nxk = nbxk(k) nyk = nbyk(k) nck = ncsk(k) npbk = npbck(k) nok = nlok(k) call dfdxk(tp,tq,npt,npk,0,nxk,nyk,nck,lxxk(1,k),lyxk(1,k), * snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k)) call dfdxk(tp(1,2),tq(1,2),npt,npk,0,nxk,nyk,nck,lxxk(1,k),lyxk(1,k), * snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k)) call zerodt(tq(1,1),nok,lok(1,k),nxk,lxxk(1,k),tp(1,3)) call zerodt(tq(1,2),nok,lok(1,k),nxk,lxxk(1,k),tp(1,3)) do j = 1, npk i = isk(j,k) abc = tq(i,3) * h(i,k) * (tq(i,1) - tmp(i,k) * tq(i,2)) pgfx(i,k) = abc fu(i,k) = fu(i,k) + abc enddo call dfdyk(tp(1,1),tq(1,1),npt,npk,0,nyk,nxk,nck,lyyk(1,k),lxyk(1,k), * snyk(1,k),isyk(1,k)) call dfdyk(tp(1,2),tq(1,2),npt,npk,0,nyk,nxk,nck,lyyk(1,k),lxyk(1,k), * snyk(1,k),isyk(1,k)) call zerodt(tq(1,1),nok,lok(1,k),nyk,lyxk(1,k),tp(1,3)) call zerodt(tq(1,2),nok,lok(1,k),nyk,lyxk(1,k),tp(1,3)) do j = 1, npk i = isk(j,k) abc = tq(i,4) * h(i,k) * (tq(i,1) - tmp(i,k) * tq(i,2) ) pgfy(i,k) = abc fv(i,k) = fv(i,k) + abc enddo if (k .lt. nz) then do j = 1, npk i = isk(j,k) dh = h(i,k) - tp(i,4) tp(i,1) = tp(i,1) + dh*(tmp(i,k) + tmp(i,k+1)) tp(i,2) = tp(i,2) + 2.*dh tp(i,4) = dh enddo endif enddo do i = 1, npt do k = 1, nzi_b(i) zfu(i) = zfu(i) + pgfx(i,k) zfv(i) = zfv(i) + pgfy(i,k) enddo enddo return end c ------------------------------------------------------------------ subroutine btpgf_old(npt,nzi_b,h,temp,dens,fu,fv,emx,emy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,lok,tp,tq,tr,lpbcwk,lpbcek,zfu,zfv, * pgfx,pgfy) c ------------------------------------------------------------------ c subroutine that calculates the pressure gradient terms in the c momentum equation for A CASE WITH CONSTANT DEPTH. c implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' include 'comm_new.h' dimension h(npt,nz),temp(npt,nz), fu(npt,nz),fv(npt,nz), * emx(1),emy(1),tp(npt,3),tq(npt,3),tr(npt,1), * dens(npt,nz),tmp(npt,nz),zfu(npt),zfv(npt), * pgfx(npt,nz),pgfy(npt,nz) dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz), * snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),nzi_b(npt), * lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz),lok(4*MAXSID,nz) pointer (p_tmp, tmp) common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) c-new temp array should contain temp, not heat content!!! c-new dens array should contain dens, not mass!!! c if (use_salt) then c .......................................case of density from EOS: c b = (grav/(1000+sigma0)) * (sigma0 - sigma(k)) c coef = -GRAVTY / (1000. + SITUD_BOT) bottom = SITUD_BOT p_tmp = loc(dens) else cc .....................case of linear (in T) density and buoyancy: c b = alph * grav * (t(k) - t_bot) coef = TALPHA * GRAVTY bottom = TEMP_BOT p_tmp = loc(temp) endif do i = 1, npt dz = h(i,1)/2. tp(i,1) = dz * tmp(i,1) tp(i,2) = dz tq(i,2) = coef * emx(i) tq(i,3) = coef * emy(i) enddo do k = 1, nz npk = nptk(k) nxk = nbxk(k) nyk = nbyk(k) nck = ncsk(k) npbk = npbck(k) nok = nlok(k) call dfdxk(tp,tq,npt,npk,0,nxk,nyk,nck,lxxk(1,k),lyxk(1,k), * snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k)) call zerodt(tq,nok,lok(1,k),nxk,lxxk(1,k),tp(1,3)) do j = 1, npk i = isk(j,k) abc = tq(i,2) * h(i,k) * tq(i,1) pgfx(i,k) = abc fu(i,k) = fu(i,k) + abc enddo call dfdyk(tp,tq,npt,npk,0,nyk,nxk,nck,lyyk(1,k),lxyk(1,k), * snyk(1,k),isyk(1,k)) call zerodt(tq,nok,lok(1,k),nyk,lyxk(1,k),tp(1,3)) do j = 1, npk i = isk(j,k) abc = tq(i,3) * h(i,k) * tq(i,1) pgfy(i,k) = abc fv(i,k) = fv(i,k) + abc enddo if (k .lt. nz) then do j = 1, npk i = isk(j,k) dz = h(i,k) - tp(i,2) tp(i,2) = dz tp(i,1) = tp(i,1) + dz * (tmp(i,k) + tmp(i,k+1)) enddo endif enddo do i = 1, npt do k = 1, nzi_b(i) zfu(i) = zfu(i) + pgfx(i,k) zfv(i) = zfv(i) + pgfy(i,k) enddo enddo return end