c: /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 c# 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(nx,ny,nz,npt,nzi,uc,vc, * w,fhd,isk,tp,h,hupi,dx2i,dy2i,ycos,inx,iny) c ------------------------------------------------------------------ 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/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) dimension isk(npt,nz),inx(nx),iny(ny) dimension uc(npt,nz),vc(npt,nz),w(npt,nz), * fhd(npt,nz), tp(npt), h(npt,nz), nzi(npt), * dx2i(1), dy2i(1), ycos(1) do i = 1, npt w(i,nz) = 0. enddo do k = 1, nz npk = nptk(k) call divergence(nx,ny,npt,npk,uc(1,k),vc(1,k),dx2i,dy2i,ycos, * isk(1,k),inx,iny,fhd(1,k)) do j = 1, npk i = isk(j,k) w(i,nz) = w(i,nz) + fhd(i,k) enddo enddo return end c ---------------------------------------------------------- subroutine wtop (npt, w, fhd, fh, h, ncyc, istep, dnt) c ---------------------------------------------------------- include 'comm_new.h' include 'diffiso.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 enddo return end c ------------------------------------------------------------------ subroutine force_u (npt,taux,tauy,fu,fv) c ------------------------------------------------------------------ dimension taux(npt), tauy(npt), fu(npt), fv(npt) do i = 1, npt fu(i) = fu(i) + taux(i) fv(i) = fv(i) + tauy(i) enddo return end c ------------------------------------------------------------------ subroutine force_u_dyice (npt,taux,tauy,cice,tauwx,tauwy,fu,fv) c ------------------------------------------------------------------ dimension taux(npt), tauy(npt), fu(npt), fv(npt) dimension tauwx(npt), tauwy(npt), cice(npt) do i = 1, npt ci = cice(i) co = 1.-ci fu(i) = fu(i) + ci* tauwx(i) +co* taux(i) fv(i) = fv(i) + ci* tauwy(i) +co* tauy(i) enddo return end c ------------------------------------------------------------------ subroutine vertu (npt,nz,nsig,nzi,bint,u,v,w,h, * fu,fv,fh) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'diffiso.h' dimension bint(1), u(npt,nz),v(npt,nz), * w(npt,nz), h(npt,nz), fu(npt,nz), fv(npt,nz), nzi(npt), fh(npt,nz) 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 w_add = w(i,k) uk = u(i,k) uk1 = u(i,k1) vk = v(i,k) vk1 = v(i,k1) ck = cupi ck1 = 1.- ck if (w_add.lt.0) then ck = ck1 ck1 = cupi endif u_ave = ck*uk + ck1*uk1 v_ave = ck*vk + ck1*vk1 u_add = h_ave*(uk1 - uk) + w_add*u_ave v_add = h_ave*(vk1 - vk) + 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 enddo enddo return end subroutine vertts (npt,nz,nzi,cint,q,qr,ep,w,h,t,ft,s,fs, * ft_va,ft_vd,fsal_va,fsal_vd,flor) c------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'diffiso.h' include 'comm_para.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) * ,ft_va(npt,nz), ft_vd(npt,nz) * ,fsal_va(npt,nz), fsal_vd(npt,nz) * ,ftp(MAXNZ),ftm(MAXNZ),fsp(MAXNZ),fsm(MAXNZ) * ,tflux(MAXNZ),sflux(MAXNZ) 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 k1 = k + 1 tk = t(i,k) tk1 = t(i,k1) sk = s(i,k) sk1 = s(i,k1) dz = 2.*h(i,k) - dz c2 = cint(k) h_ave = c2/dz t_dif = h_ave*(tk1 - tk) s_dif = h_ave*(sk1 - sk) ft_vd(i,k) = ft_vd(i,k) + flor*t_dif ft_vd(i,k1) = ft_vd(i,k1) - flor*t_dif fsal_vd(i,k) = fsal_vd(i,k) + flor*s_dif fsal_vd(i,k1) = fsal_vd(i,k1) - flor*s_dif ft(i,k) = ft(i,k) + t_dif ft(i,k1) = ft(i,k1) - t_dif fs(i,k) = fs(i,k) + s_dif fs(i,k1) = fs(i,k1) - s_dif enddo enddo if (use_tvd_vert) then do i = 1, npt nzc = nzi(i) ftp(nzc) = 0. ftm( 1) = 0. fsp(nzc) = 0. fsm( 1) = 0. do k = 1, nzc - 1 wik = w(i,k) kp = k + 1 tk = t(i,k) tkp = t(i,kp) sk = s(i,k) skp = s(i,kp) wp = 0. if (wik.gt.0.)wp = wik wm = 0. if (wik.lt.0.)wm = wik tflux(k) = wm* tk + wp* tkp sflux(k) = wm* sk + wp* skp dt = tk - tkp ds = sk - skp ftp(k ) = wp*dt ftm(kp) = wm*dt fsp(k ) = wp*ds fsm(kp) = wm*ds enddo do k = 1, nzc - 1 kp = k + 1 rp = 0 if (ftp(kp).ne.0) rp = ftp(k)/ftp(kp) rp = max(0.,rp) rm = 0 if (ftm(k ).ne.0) rm = ftm(kp)/ftm(k) rm = max(0.,rm) sp = rp/(1.+rp) sm = rm/(1.+rm) tflux(k) = tflux(k) + sp*ftp(kp)- sm*ftm(k) rp = 0 if (fsp(kp).ne.0) rp = fsp(k)/fsp(kp) rp = max(0.,rp) rm = 0 if (fsm(k ).ne.0) rm = fsm(kp)/fsm(k) rm = max(0.,rm) sp = rp/(1.+rp) sm = rm/(1.+rm) sflux(k) = sflux(k) + sp*fsp(kp)- sm*fsm(k) t_adv = tflux(k) s_adv = sflux(k) ft_va(i,k) = ft_va(i,k) + flor*t_adv ft_va(i,kp) = ft_va(i,kp) - flor*t_adv fsal_va(i,k) = fsal_va(i,k) + flor*s_adv fsal_va(i,kp) = fsal_va(i,kp) - flor*s_adv ft(i,k) = ft(i,k) + t_adv ft(i,kp) = ft(i,kp) - t_adv fs(i,k) = fs(i,k) + s_adv fs(i,kp) = fs(i,kp) - s_adv enddo enddo else do i = 1, npt do k = 1, nzi(i) - 1 k1 = k + 1 tk = t(i,k) tk1 = t(i,k1) sk = s(i,k) sk1 = s(i,k1) wik = w(i,k) ck = cupi_ts ck1 = 1.- ck if (wik.lt.0) then ck = ck1 ck1 = cupi_ts endif t_ave = ck*tk + ck1*tk1 s_ave = ck*sk + ck1*sk1 t_adv = wik*t_ave s_adv = wik*s_ave ft_va(i,k) = ft_va(i,k) + flor*t_adv ft_va(i,k1) = ft_va(i,k1) - flor*t_adv fsal_va(i,k) = fsal_va(i,k) + flor*s_adv fsal_va(i,k1) = fsal_va(i,k1) - flor*s_adv ft(i,k) = ft(i,k) + t_adv ft(i,k1) = ft(i,k1) - t_adv fs(i,k) = fs(i,k) + s_adv fs(i,k1) = fs(i,k1) - s_adv enddo enddo endif 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,ntr,nzi,cint,w,h,tr,ftr,ftr_va,ftr_vd,flor) c------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'diffiso.h' include 'comm_para.h' dimension cint(nz), w(npt,nz), h(npt,nz), tr(npt,nz,1), ftr(npt,nz,1) dimension ftr_va(npt,nz+1,ntr), ftr_vd(npt,nz+1,ntr) dimension nzi(npt) parameter (H_ATT1 = -1./17.) dimension ftp(MAXNZ),ftm(MAXNZ),tflux(MAXNZ) if (use_tvd_vert) then do n = 1, ntr do i = 1, npt nzc = nzi(i) ftp(nzc) = 0. ftm( 1) = 0. do k = 1, nzc - 1 wik = w(i,k) kp = k + 1 tk = tr(i,k,n) tkp = tr(i,kp,n) wp = 0. if (wik.gt.0.)wp = wik wm = 0. if (wik.lt.0.)wm = wik tflux(k) = wm* tk + wp* tkp dt = tk - tkp ftp(k ) = wp*dt ftm(kp) = wm*dt enddo do k = 1, nzc - 1 kp = k + 1 rp = 0 if (ftp(kp).ne.0) rp = ftp(k)/ftp(kp) rp = max(0.,rp) rm = 0 if (ftm(k ).ne.0) rm = ftm(kp)/ftm(k) rm = max(0.,rm) sp = rp/(1.+rp) sm = rm/(1.+rm) tflux(k) = tflux(k) + sp*ftp(kp)- sm*ftm(k) t_adv = tflux(k) ftr_va(i,k,n) = ftr_va(i,k,n) + flor*t_adv ftr_va(i,kp,n) = ftr_va(i,kp,n) + flor*t_adv ftr(i,k,n) = ftr(i,k,n) + t_adv ftr(i,kp,n) = ftr(i,kp,n) - t_adv enddo enddo enddo else do i = 1, npt do k = 1, nzi(i) - 1 k1 = k + 1 wik = w(i,k) ck = cupi_tr ck1 = 1.- ck if (wik.lt.0) then ck = ck1 ck1 = cupi_tr endif do n = 1, ntr trk = tr(i,k,n) trk1 = tr(i,k1,n) tr_ave = ck*trk + ck1*trk1 tr_adv = if_tr_adv*wik*tr_ave ftr_va(i,k,n) = ftr_va(i,k,n) + flor*tr_adv ftr_va(i,k1,n) = ftr_va(i,k1,n) - flor*tr_adv ftr(i,k,n) = ftr(i,k,n) + tr_adv ftr(i,k1,n) = ftr(i,k1,n) - tr_adv enddo enddo enddo endif do i = 1, npt dz = h(i,1) do k = 1, nzi(i) - 1 k1 = k + 1 hval = h(i,k) c2 = cint(k) dz = 2.*hval - dz h_ave = c2/dz do n = 1, ntr trk = tr(i,k,n) trk1 = tr(i,k1,n) tr_dif = h_ave*(trk1 - trk) ftr_vd(i,k,n) = ftr_vd(i,k,n) + flor*tr_dif ftr_vd(i,k1,n) = ftr_vd(i,k1,n) - flor*tr_dif ftr(i,k,n) = ftr(i,k,n) + tr_dif ftr(i,k1,n) = ftr(i,k1,n) - tr_dif 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(npt,nz),en(npten,1), * nptk(1),isk(npt,1),basin(nz) 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,k) * h(i,k) * (uu*uu + vv*vv) enddo en(1,k)= 0.5 * ek / basin(k) 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(npt,1),t(npt,1),dens(npt,1),w(npt,1),basin(nz), * en(npten,1), tp(1), tmp(npt,1), nptk(1), isk(npt,1) pointer (ptmp, tmp) ptmp = loc(dens) base = SITUD_BOT shift = 1000. 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,k) * 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(t_fac,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)/t_fac vc(i,k) = vc(i,k) + binv*fv(i,k)/t_fac fu(i,k) = abinv*fu(i,k) fv(i,k) = abinv*fv(i,k) enddo enddo return end c ------------------------------------------------------------------ subroutine bc_steer(npt,fu,fv,nzi,lxxk,lyyk,iox,mask,isteer,saxpk,saypk) 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 iox(npt), mask(nxp*nyp,nz) dimension lxxk(MXBDY,nz), lyyk(MXBDY,nz), saxpk(npt,nz), saypk(npt,nz) dimension nzi(npt) dimension fu(npt,nz),fv(npt,nz) c redistribute baroclinic momentum forcings to "steer the flow" c near boundaries if (isteer.eq.1) then nxy = nxp*nyp do k = 1, nz do i = 1, nbxk(k) ik = lxxk(i,k) gridfac = saypk(ik,k)/saxpk(ik,k) ixy = iox(ik) i_N = 0 i_S = 0 if (ixy+nxp.le.nxy) i_N = mask(ixy + nxp,k) if (ixy-nxp.ge.1) i_S = mask(ixy - nxp,k) if (i_N.ne.0) fv(i_N,k) = fv(i_N,k) - 0.5 * fu(ik,k)* gridfac if (i_S.ne.0) fv(i_S,k) = fv(i_S,k) + 0.5 * fu(ik,k)* gridfac fu(ik,k) = 0. enddo do i = 1, nbyk(k) ik = lyyk(i,k) gridfac = saxpk(ik,k)/saypk(ik,k) ixy = iox(ik) i_E = 0 i_W = 0 if (ixy+1.le.nxy) i_E = mask(ixy + 1,k) if (ixy-1.ge.1) i_W = mask(ixy - 1,k) if (i_E.ne.0) fu(i_E,k) = fu(i_E,k) - 0.5 * fv(ik,k)* gridfac if (i_W.ne.0) fu(i_W,k) = fu(i_W,k) + 0.5 * fv(ik,k)* gridfac fv(ik,k) = 0. enddo enddo endif return c end of bc_relax end c ------------------------------------------------------------------ subroutine horizu(nx,ny,nz,npt,nsig,u,v,uc,vc,f,fu,fv,fhd,tp,isk, * lxxk,lyyk,nzibx,nziby,corx,cory,fh,hupi,dx2i,dy2i,ycos,inx,iny) 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) transport arrays c isk = (input) indices to convert from compressed to full c tp = (input) temporary space. implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' common/friction/ ibfric, b_fric, sw_fric common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) dimension isk(npt,nz),inx(nx),iny(ny),lxxk(MXBDY,nz),lyyk(MXBDY,nz) dimension u(npt,nz),v(npt,nz),uc(npt,nz),vc(npt,nz), * f(npt),fu(npt,nz),fv(npt,nz), * tp(npt),fh(npt,nz),nzibx(npt),nziby(npt),fhd(npt,nz), * corx(npt,nz), cory(npt,nz), * dx2i(1), dy2i(1), ycos(1) 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 c.....compute the horizontal terms of the momentum forcings, fu,fv do k = 1, nz npk = nptk(k) call divergence_up(nx,ny,npt,npk,u(1,k),uc(1,k),vc(1,k), * dx2i,dy2i,ycos,isk(1,k),inx,iny,tp,hupi) do j = 1, npk i = isk(j,k) c.....add coriolis terms: tmpu = f(i)*vc(i,k) fu(i,k) = fu(i,k) + tmpu corx(i,k) = tmpu tmpv = -f(i)*uc(i,k) fv(i,k) = fv(i,k) + tmpv cory(i,k) = tmpv c.....add momentum advection terms: fu(i,k) = fu(i,k) - tp(i) enddo call divergence_up(nx,ny,npt,npk,v(1,k),uc(1,k),vc(1,k), * dx2i,dy2i,ycos,isk(1,k),inx,iny,tp,hupi) do j = 1, npk i = isk(j,k) c.....add momentum advection terms: fv(i,k) = fv(i,k) - tp(i) enddo enddo if (ibfric.ne.0) then c..... add bottom drag in last layer, and along horizontal boundaries do i = 1, npt mz = nzibx(i) if (mz.ne.0) then bfricu = b_fric * u(i,mz) if (ibfric.eq.2) then bfricu = bfricu* abs(u(i,mz)) endif fu(i,mz) = fu(i,mz) - bfricu endif mz = nziby(i) if (mz.ne.0) then bfricv = b_fric * v(i,mz) if (ibfric.eq.2) then bfricv = bfricv* abs(v(i,mz)) endif fv(i,mz) = fv(i,mz) - bfricv endif enddo if (sw_fric.ne.0) then do k = 1, nz do ib = 1, nbxk(k) i = lxxk(ib,k) bfricv = sw_fric * v(i,k) if (ibfric.eq.2) then bfricv = bfricv* abs(v(i,k)) endif fv(i,k) = fv(i,k) - bfricv enddo do ib = 1, nbyk(k) i = lyyk(ib,k) bfricu = sw_fric * u(i,k) if (ibfric.eq.2) then bfricu = bfricu* abs(u(i,k)) endif fu(i,k) = fu(i,k) - bfricu enddo enddo endif endif return c end of horizu. end c ------------------------------------------------------------------ subroutine horizt(nx,ny,nz,npt,h,uc,vc,u,v,t,ft,isk,tp,if_adv, * ft_ha,flor,tupi,dx2i,dy2i,ycos,inx,iny) c ------------------------------------------------------------------ c implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' include 'diffiso.h' 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), * u(npt,nz),v(npt,nz), tp(npt),h(npt,nz),ft_ha(npt,nz), * dx2i(1), dy2i(1), ycos(1) dimension isk(npt,nz),inx(nx),iny(ny) if (if_adv.eq.0) return do k = 1, nz npk = nptk(k) if (use_tvd_horz) then call divergence_tvd(nx,ny,npt,npk,t(1,k),uc(1,k),vc(1,k), * dx2i,dy2i,ycos,isk(1,k),inx,iny,tp,tupi) else call divergence_up(nx,ny,npt,npk,t(1,k),uc(1,k),vc(1,k), * dx2i,dy2i,ycos,isk(1,k),inx,iny,tp,tupi) endif do j = 1, npk i = isk(j,k) f_adv = tp(i) ft(i,k) = ft(i,k) - f_adv ft_ha(i,k) = ft_ha(i,k) - flor*f_adv enddo enddo return c end of horizt 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 tupdat(t_fac,npt,nz,nzi,binv,abinv,h,t,ft) c ------------------------------------------------------------------ c update temperature fields as was done for u, v in vel_updat c implicit real(a-h,o-z),integer(i-n) dimension t(npt,nz),ft(npt,nz),nzi(npt),h(npt,nz),t_fac(nz) c do i = 1, npt do k = 1, nzi(i) t(i,k) = t(i,k) + binv*ft(i,k)/h(i,k)/t_fac(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) 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) 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 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 ------------------------------------------------------------------ subroutine bc_relax(dt,coef,mbc,lxxk,lyyk,npt,u,v,fu,fv,nzi) c ------------------------------------------------------------------ c relax to 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) dimension fu(npt,nz),fv(npt,nz) 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 const = coef/dt if(mbc.eq.1 .or. mbc.eq.4) then do k=1,nz do i=1,nbyk(k)+ncsk(k) ik = lyyk(i,k) fu(ik,k) = fu(ik,k) - const* u(ik,k) enddo enddo endif c if(mbc.eq.1 .or. mbc.eq.3) then do k=1,nz do i=1,nbxk(k)+ncsk(k) ik = lxxk(i,k) fv(ik,k) = fv(ik,k) - const* v(ik,k) enddo enddo endif return c end of bc_relax end c -------------------------------------------------------- logical function non_stable (iout, npt, nxp, nz, nzi, 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), nzi(npt) 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 do k = 1, nzi(i) if (t(i,k) .lt. -10..or.t(i,k).gt.50) then jj = 1 + (iox(i)-1)/nxp ii = iox(i) - (jj-1)*nxp write (iout, 11) i, ii, jj, k, t(i,2) endif if (u(i,k)**2 + v(i,k)**2 .gt. 400.) then jj = 1 + (iox(i)-1)/nxp ii = iox(i) - (jj-1)*nxp write (iout, 12) i, ii, jj, k endif enddo enddo endif 11 format ('t(k=',i4,',2)[i=',i3,',j=',i3,',k=',i3,'] =', g10.3) 12 format ('i,ii,jj,k=',4i8,' 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 pgfu(npt,nzi,h,dens,fu,fv,saxmk,saxpk,saymk,saypk, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,lok,tp,dpres,ddens, * lpbcwk,lpbcek,pgfx,pgfy,gradE) 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' real*8 dpres(npt),ddens(npt) dimension h(npt,nz),fu(npt,nz),fv(npt,nz), * tp(npt,4), * saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz), * dens(npt,nz), * pgfx(npt,nz),pgfy(npt,nz) dimension gradE(npt,2) 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(npt), * lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz),lok(4*MAXSID,nz) 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 dens array should contain dens, not mass!!! c c .......................................case of density from EOS: c b = (grav/(1000+sigma0)) * (sigma0 - sigma(k)) coef = -GRAVTY/ (1000. + SITUD_BOT) do i = 1, npt dh = h(i,1)/2. dpres(i) = dh*dens(i,1) ddens(i)= dens(i,1) tp(i,2) = dh tp(i,4) = dh enddo do k = 1, nz kp = k+1 npk = nptk(k) nxk = nbxk(k) nyk = nbyk(k) nck = ncsk(k) npbk = npbck(k) nok = nlok(k) call dfdx_dp(dpres,tp(1,1),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), * saxmk(1,k),saxpk(1,k)) call dfdx_dp(ddens,tp(1,3),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), * saxmk(1,k),saxpk(1,k)) do j = 1, npk i = isk(j,k) abc = coef * h(i,k) * tp(i,1) pgfx(i,k) = abc fu(i,k) = fu(i,k) + abc gradE(i,1) = gradE(i,1) + coef*h(i,k)*tp(i,2)*tp(i,3) enddo call dfdy_dp(dpres,tp(1,1),npt,npk,0,nyk,nxk,nck,lyyk(1,k), * lxyk(1,k),snyk(1,k),isyk(1,k), * saymk(1,k),saypk(1,k)) call dfdy_dp(ddens,tp(1,3),npt,npk,0,nyk,nxk,nck,lyyk(1,k), * lxyk(1,k),snyk(1,k),isyk(1,k), * saymk(1,k),saypk(1,k)) do j = 1, npk i = isk(j,k) abc = coef * h(i,k) * tp(i,1) pgfy(i,k) = abc fv(i,k) = fv(i,k) + abc gradE(i,2) = gradE(i,2) + coef*h(i,k)*tp(i,2)*tp(i,3) enddo if (k .lt. nz) then do j = 1, npk i = isk(j,k) dh = h(i,k) - tp(i,4) dpres(i) = dpres(i) + dh*(dens(i,k) + dens(i,kp)) ddens(i)= dens(i,kp) tp(i,2) = tp(i,2) + 2.*dh tp(i,4) = dh enddo endif enddo return end c ------------------------------------- subroutine tdecap(npt, nz, nsig, t, h) c ------------------------------------- c.....convert heat content to temperature. dimension t(npt,nz), h(npt,nz) do i = 1, npt do k = 1, nsig t(i,k) = t(i,k)/h(i,k) enddo enddo return end c ------------------------------------------------------------------ subroutine capt(npt, nz, nsig, t, h) c ------------------------------------------------------------------ c form the heat content, ht, from the temperature. c dimension t(npt,nz),h(npt,nz) c do i = 1, npt do k = 1, nsig t(i,k) = t(i,k)*h(i,k) enddo enddo return c end of capt. end c ------------------------------------------------------------------ subroutine flux_ca(nptz,dtmix,t,h,t_old,ft_ca) c ------------------------------------------------------------------ dimension t(nptz),h(nptz),t_old(nptz),ft_ca(nptz) do i = 1, nptz ft_ca(i) = ft_ca(i) + h(i) * ( t(i) - t_old(i) ) / dtmix enddo return end c ------------------------------------------------------------------ subroutine bndrys_new (npt,nzi,lxxk,lyyk,nzibx,nziby) c ------------------------------------------------------------------ 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 nzi(npt), nzibx(npt), nziby(npt) do i = 1, npt nzibx(i) = nzi(i) nziby(i) = nzi(i) enddo do k = nz, 1, -1 do ib = 1, nbxk(k) i = lxxk(ib,k) nzibx(i) = nzibx(i) - 1 enddo do ib = 1, nbyk(k) i = lyyk(ib,k) nziby(i) = nziby(i) - 1 enddo enddo return end