subroutine HDIFSN (hfu, hfv, hft, hfr) !----------------------------------------------------------------------- ! ! HDIFSN computes horizontal diffusion tendency using a non-linear ! smagorinsky scheme. the final result are diffusive flux ! divergences for u, v, t and r (hfu, hfv, hft and hfr). these ! divergences will be added to vertical diffusion to obtain total ! diffusion tendency. ! !----------------------------------------------------------------------- use allocate_mod include 'RESOLUTION.h' include 'PARAMETERS.h' include 'STDUNITS.h' include 'BNFX.h' include 'CLCLE.h' include 'CONMLEV.h' include 'CONSLEV.h' include 'FILEC.h' include 'FILEPN.h' include 'FILEPC.h' include 'FILEPS.h' include 'GDINFO.h' include 'HDICUM.h' include 'LIMIT.h' include 'NJJJ.h' include 'SET.h' !----------------------------------------------------------------------- ! user interface variables !----------------------------------------------------------------------- real, intent (out), dimension(kmax,ibeg:iend) :: hfu real, intent (out), dimension(kmax,ibeg:iend) :: hfv real, intent (out), dimension(kmax,ibeg:iend) :: hft real, intent (out), dimension(kmax,ibeg:iend) :: hfr !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- real, dimension(imx) :: coxw real, dimension(imx) :: coyn1 real, dimension(imx) :: coys1 real, dimension(imx) :: csrln real, dimension(imx) :: csrls real, dimension(imx) :: gaman real, dimension(imx) :: gamas real, dimension(imx) :: psci real, dimension(imx) :: pscni real, dimension(imx) :: pscsi real, dimension(imx) :: psni real, dimension(imx) :: pssi real, dimension(imx) :: sqk0c real, dimension(imx) :: sqk0n real, dimension(imx) :: sqk0s real, dimension(imx) :: sqkci real, dimension(imx) :: sqkcsi real, dimension(imx) :: sqkni real, dimension(imx) :: sqksi real, dimension(imx+1) :: gamaw real, dimension(imx+1) :: psbarw real, dimension(imx+1) :: rlxw real, dimension(imx+1) :: sqkw real, dimension(kmax,imx) :: csrltn real, dimension(kmax,imx) :: csrlts real, dimension(kmax,imx) :: dfn1 real, dimension(kmax,imx) :: dfs1 real, dimension(kmax,imx) :: drn1 real, dimension(kmax,imx) :: drs1 real, dimension(kmax,imx) :: dtn1 real, dimension(kmax,imx) :: dts1 real, dimension(kmax,imx) :: dun1 real, dimension(kmax,imx) :: dus1 real, dimension(kmax,imx) :: dvn1 real, dimension(kmax,imx) :: dvs1 real, dimension(kmax,imx) :: hryn real, dimension(kmax,imx) :: hryn1 real, dimension(kmax,imx) :: hrys real, dimension(kmax,imx) :: hrys1 real, dimension(kmax,imx) :: htyn real, dimension(kmax,imx) :: htyn1 real, dimension(kmax,imx) :: htys real, dimension(kmax,imx) :: htys1 real, dimension(kmax,imx) :: rci real, dimension(kmax,imx) :: rcni real, dimension(kmax,imx) :: rcsi real, dimension(kmax,imx) :: reyn1 real, dimension(kmax,imx) :: reys1 real, dimension(kmax,imx) :: rlynn real, dimension(kmax,imx) :: rlyss real, dimension(kmax,imx) :: rni real, dimension(kmax,imx) :: rsi real, dimension(kmax,imx) :: sxyn real, dimension(kmax,imx) :: sxyn1 real, dimension(kmax,imx) :: sxys real, dimension(kmax,imx) :: sxys1 real, dimension(kmax,imx) :: syyn real, dimension(kmax,imx) :: syyn1 real, dimension(kmax,imx) :: syys real, dimension(kmax,imx) :: syys1 real, dimension(kmax,imx) :: tci real, dimension(kmax,imx) :: tcni real, dimension(kmax,imx) :: tcsi real, dimension(kmax,imx) :: tni real, dimension(kmax,imx) :: tsi real, dimension(kmax,imx) :: uci real, dimension(kmax,imx) :: ucni real, dimension(kmax,imx) :: ucsi real, dimension(kmax,imx) :: uni real, dimension(kmax,imx) :: usi real, dimension(kmax,imx) :: vci real, dimension(kmax,imx) :: vcni real, dimension(kmax,imx) :: vcsi real, dimension(kmax,imx) :: vni real, dimension(kmax,imx) :: vsi real, dimension(kmax,imx) :: wtn real, dimension(kmax,imx) :: wts real, dimension(kmax,imx+1) :: dfw real, dimension(kmax,imx+1) :: drw real, dimension(kmax,imx+1) :: dtw real, dimension(kmax,imx+1) :: duw real, dimension(kmax,imx+1) :: dvw real, dimension(kmax,imx+1) :: hrxw real, dimension(kmax,imx+1) :: htxw real, dimension(kmax,imx+1) :: sxxw real, dimension(kmax,imx+1) :: syxw real, dimension(kmax,imx+1) :: wtw !----------------------------------------------------------------------- ! internal variables sith storage association !----------------------------------------------------------------------- save cumdy real, dimension (2,nfr) :: cumdy !----------------------------------------------------------------------- ! code !----------------------------------------------------------------------- j = IFIX(tjlopc(i1)) do i = 1,imx+1 gamaw(i) = .5 enddo if (nstflg .EQ. 1) then do i = i1,i2 + 1 rlxw(i) = 1.0e - 9 enddo do i = i1,i2+ngr do k = 1,kmax sxxw(k,i) = 0.0 syxw(k,i) = 0.0 htxw(k,i) = 0.0 hrxw(k,i) = 0.0 enddo enddo endif if ((j .LE. 2 .OR. j .GE. jmax - 1) .AND. nstflg .EQ. 1) then do i = i1,i2 do k = 1,kmax sxys (k,i) = 0.0 syys (k,i) = 0.0 htys (k,i) = 0.0 hrys (k,i) = 0.0 sxyn (k,i) = 0.0 syyn (k,i) = 0.0 htyn (k,i) = 0.0 hryn (k,i) = 0.0 enddo enddo do i = i1,i2+ngr do k = 1,kmax sxys1(k,i) = 0.0 syys1(k,i) = 0.0 htys1(k,i) = 0.0 hrys1(k,i) = 0.0 sxyn1(k,i) = 0.0 syyn1(k,i) = 0.0 htyn1(k,i) = 0.0 hryn1(k,i) = 0.0 enddo enddo endif do i = 1,imax sqk0c(i) = .04 sqk0s(i) = .04 sqk0n(i) = .04 enddo if (nstflg .EQ. 1 .AND. i1 .EQ. 1.) then do k = 1,kmax !----------------------------------------------------------------------- ! sxxw(k,1) = bnfxm2(k,1,jj,21, ngd-1) ! syxw(k,1) = bnfxm2(k,1,jj,22, ngd-1) ! htxw(k,1) = bnfxm2(k,1,jj,23,ngd-1) ! hrxw(k,1) = bnfxm2(k,1,jj,24,ngd-1) !----------------------------------------------------------------------- sxxw(k,1) = bnfxm2(k,1,jj,27, ngd-1) syxw(k,1) = bnfxm2(k,1,jj,28, ngd-1) htxw(k,1) = bnfxm2(k,1,jj,29,ngd-1) hrxw(k,1) = bnfxm2(k,1,jj,30,ngd-1) enddo endif if (nstflg .EQ. 1 .AND. i2 .EQ. imax) then do k = 1,kmax !----------------------------------------------------------------------- ! sxxw(k,imax + 1) = bnfxm2(k,2,jj,21, ngd-1) ! syxw(k,imax + 1) = bnfxm2(k,2,jj,22, ngd-1) ! htxw(k,imax + 1) = bnfxm2(k,2,jj,23,ngd-1) ! hrxw(k,imax + 1) = bnfxm2(k,2,jj,24,ngd-1) !----------------------------------------------------------------------- sxxw(k,imax + 1) = bnfxm2(k,2,jj,27, ngd-1) syxw(k,imax + 1) = bnfxm2(k,2,jj,28, ngd-1) htxw(k,imax + 1) = bnfxm2(k,2,jj,29,ngd-1) hrxw(k,imax + 1) = bnfxm2(k,2,jj,30,ngd-1) enddo endif !----------------------------------------------------------------------- ! initialize arrays for storing boundary fluxes at mesh interface ! inner edge of window frame !----------------------------------------------------------------------- if (nstflg .EQ. 0) go to 100 if (jrow .NE. 1) go to 100 if (i1 .EQ. 3) then cumdy(1,jj) = 0. do iid = 1 , 4 do k = 1,kmax cumhid(k,1,jj,iid) = 0. enddo enddo endif if (i2 .EQ. imax - 2) then cumdy(2,jj) = 0. do iid = 1 , 4 do k = 1,kmax cumhid(k,2,jj,iid) = 0. enddo enddo endif 100 continue incc = 1 incn = 1 incs = 1 if (nstflg .EQ. 0) go to 160 !----------------------------------------------------------------------- ! nested case !----------------------------------------------------------------------- if (i1 .GT. 2 .AND. i1 .LT. imax - 1) then if (j .LE. 2 .OR. j .GE. jmax - 1) incc = ngr if (j .LE. 1 .OR. j .GE. jmax - 2) incn = ngr if (j .LE. 3 .OR. j .GE. jmax ) incs = ngr endif !----------------------------------------------------------------------- ! interpolate to get fluxes at east and west mesh interfaces of ! central section !----------------------------------------------------------------------- if (jrow .EQ. njrow) then if (i2 .EQ. 2) then do k = 1,kmax sxxw(k,3) = cumhid(k,1,jj,1) syxw(k,3) = cumhid(k,1,jj,2) htxw(k,3) = cumhid(k,1,jj,3) hrxw(k,3) = cumhid(k,1,jj,4) enddo endif if (i1 .EQ. imax - 1) then i = imax - 1 do k = 1,kmax sxxw(k,i) = cumhid(k,2,jj,1) syxw(k,i) = cumhid(k,2,jj,2) htxw(k,i) = cumhid(k,2,jj,3) hrxw(k,i) = cumhid(k,2,jj,4) enddo endif endif if (i1 .EQ. 3) then gamaw(i1) = dlonpc(i1)/(dlonpc(i1) + dlonpc(i1 - 1)) rlxw (i1) = 1./(.5*(dlonpc(i1) + dlonpc(i1 - 1))* & arad*COS(rlatpc(i1))) if (rwdrt .GE. .5) go to 120 d1 = .5 - rwdrt d2 = 1. - d1 psbarw(i1) = (d1*psps (i1 - 1) + d2*pspc (i1 - 1))*gamaw(i1) + & pspc (i1)*(1. - gamaw(i1)) sqkw(i1) = (d1*sqk0s(i1 - 1) + d2*sqk0c(i1 - 1))*gamaw(i1) + & sqk0c(i1)*(1. - gamaw(i1)) do k = 1,kmax duw(k,i1) = (upc(k,i1) - (d1*ups(k,i1 - 1) + & d2*upc(k,i1 - 1)))*rlxw(i1) - TAN(rlatpc(i1))* & ((d1*vps(k,i1 - 1) + d2*vpc(k,i1 - 1))* & gamaw(i1) + vpc(k,i1)*(1. - gamaw(i1)))/arad dvw(k,i1) = (vpc(k,i1) - (d1*vps(k,i1 - 1) + & d2*vpc(k,i1 - 1)))*rlxw(i1) dtw(k,i1) = (tpc(k,i1) - (d1*tps(k,i1 - 1) + & d2*tpc(k,i1 - 1)))*rlxw(i1) drw(k,i1) = (rpc(k,i1) - (d1*rps(k,i1 - 1) + & d2*rpc(k,i1 - 1)))*rlxw(i1) enddo go to 140 120 d1 = 1.5 - rwdrt d2 = 1. - d1 psbarw(i1) = (d1*pspc (i1 - 1) + d2*pspn (i1 - 1))*gamaw(i1) + & pspc (i1)*(1. - gamaw(i1)) sqkw(i1) = (d1*sqk0c(i1 - 1) + d2*sqk0n(i1 - 1))*gamaw(i1) + & sqk0c(i1)*(1. - gamaw(i1)) do k = 1,kmax duw(k,i1) = (upc(k,i1) - (d1*upc(k,i1 - 1) + & d2*upn(k,i1 - 1)))*rlxw(i1) - TAN(rlatpc(i1))* & ((d1*vpc(k,i1 - 1) +d2*vpn(k,i1 - 1))* & gamaw(i1) + vpc(k,i1)*(1. - gamaw(i1)))/arad dvw(k,i1) = (vpc(k,i1) - (d1*vpc(k,i1 - 1) + & d2*vpn(k,i1 - 1)))*rlxw(i1) dtw(k,i1) = (tpc(k,i1) - (d1*tpc(k,i1 - 1) + & d2*tpn(k,i1 - 1)))*rlxw(i1) drw(k,i1) = (rpc(k,i1) - (d1*rpc(k,i1 - 1) + & d2*rpn(k,i1 - 1)))*rlxw(i1) enddo endif 140 if (i2 .EQ. imax - 2) then i2p = i2 + 1 i2m = i2p - incc gamaw(i2p) = dlonpc(i2p)/(dlonpc(i2p) + dlonpc(i2m)) rlxw (i2p) = 1./(.5*(dlonpc(i2p) + dlonpc(i2m))* & arad*COS(rlatpc(i2m))) if (rwdrt .GE. .5) go to 150 d1 = .5 - rwdrt d2 = 1. - d1 psbarw(i2p) = (d1*psps(i2p) + d2*pspc(i2p))* & (1. - gamaw(i2p)) + pspc (i2m)*gamaw(i2p) sqkw(i2p) = (d1*sqk0s(i2p) + d2*sqk0c(i2p))* & (1. - gamaw(i2p)) + sqk0c(i2m)*gamaw(i2p) do k = 1,kmax duw(k,i2p) = ((d1*ups(k,i2p) + d2*upc(k,i2p)) - upc(k,i2m))* & rlxw(i2p) - TAN(rlatpc(i2m))*((d1*vps(k,i2p) + & d2*vpc(k,i2p))*(1. - gamaw(i2p)) + vpc(k,i2m)* & gamaw(i2p))/arad dvw(k,i2p) = ((d1*vps(k,i2p) + d2*vpc(k,i2p)) - vpc(k,i2m))* & rlxw(i2p) dtw(k,i2p) = ((d1*tps(k,i2p) + d2*tpc(k,i2p)) - tpc(k,i2m))* & rlxw(i2p) drw(k,i2p) = ((d1*rps(k,i2p) + d2*rpc(k,i2p)) - rpc(k,i2m))* & rlxw(i2p) enddo go to 160 150 d1 = 1.5 - rwdrt d2 = 1. - d1 psbarw(i2p) = (d1*pspc(i2p) + d2*pspn(i2p))* & (1. - gamaw(i2p)) + pspc(i2m)*gamaw(i2p) sqkw(i2p) = (d1*sqk0c(i2p) + d2*sqk0n(i2p))*(1. - gamaw(i2p)) + & sqk0c(i2m)*gamaw(i2p) do k = 1,kmax duw(k,i2p) = ((d1*upc(k,i2p) + d2*upn(k,i2p)) - upc(k,i2m))* & rlxw(i2p) - TAN(rlatpc(i2m))*((d1*vpc(k,i2p) + & d2*vpn(k,i2p))*(1. - gamaw(i2p)) + vpc(k,i2m)* & gamaw(i2p))/arad dvw(k,i2p) = ((d1*vpc(k,i2p) + d2*vpn(k,i2p)) - vpc(k,i2m))* & rlxw(i2p) dtw(k,i2p) = ((d1*tpc(k,i2p) + d2*tpn(k,i2p)) - tpc(k,i2m))* & rlxw(i2p) drw(k,i2p) = ((d1*rpc(k,i2p) + d2*rpn(k,i2p)) - rpc(k,i2m))* & rlxw(i2p) enddo endif 160 continue !----------------------------------------------------------------------- ! outer grid begins here !----------------------------------------------------------------------- do i = i1,i2 stpc(i) = stpc(i) + deltat enddo incr = MIN0(incn,incc) inc1 = MAX0(incc,incn) inc2 = (inc1 + 1)/2 insc = MIN0(incs,incc) insc1 = MAX0(incc,incs) insc2 = (insc1 + 1)/2 !----------------------------------------------------------------------- ! compute weights for interpolation to interface "gaman, rlynn" ! and weights for interface fluxes "wtw, wtn, wts" !----------------------------------------------------------------------- csrlt = COS(rlatpc(i1)) csrlnn = COS(rlatpn(i1)) csrlss = COS(rlatps(i1)) tnrlt = TAN(rlatpc(i1)) wtdty = arad*dlatpc(i1)/dlspc(i1) do i = i1,i2 csrls(i) = COS(rlatpc(i) - dlatpc(i)*.5) csrln(i) = COS(rlatpc(i) + dlatpc(i)*.5) enddo do i = i1,i2 do k = 1,kmax csrlts(k,i) = csrls(i) csrltn(k,i) = csrln(i) enddo enddo wtn1 = (arad*COS(rlatpc(i1) + dlatpc(i1)*.5)*dlonpc(i1))/ & dlspc(i1) wts1 = (arad*COS(rlatpc(i1) - dlatpc(i1)*.5)*dlonpc(i1))/ & dlspc(i1) gaman1 = dlatpn(i1)/(dlatpn(i1) + dlatpc(i1)) gamas1 = dlatps(i1)/(dlatps(i1) + dlatpc(i1)) rlyn1 = 1./(.5*(dlatpn(i1) + dlatpc(i1))*arad) rlys1 = 1./(.5*(dlatps(i1) + dlatpc(i1))*arad) do i = i1,i2 do k = 1,kmax rlynn(k,i) = rlyn1 rlyss(k,i) = rlys1 wtw(k,i) = wtdty wtn(k,i) = wtn1 wts(k,i) = wts1 enddo enddo do i = i1,i2 gaman(i) = gaman1 gamas(i) = gamas1 enddo !----------------------------------------------------------------------- ! if central and north rows are different resolution (incc.ne.incn) ! then interpolate coarser row to finer resolution, storing result ! in "xni,xci; x=u,v,t,r" !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! fluxes for inner nest for j=2 or j=jmax-2 !----------------------------------------------------------------------- if (incc .NE. incn) then do ii = 1,inc1 dist = FLOAT(2*ii - 1)/FLOAT(2*inc1) if (ii .GT. inc2) go to 220 ip = 0 im = inc1 d1 = .5 + dist d2 = .5 - dist go to 230 220 continue ip = inc1 im = 0 d1 = - .5 + dist d2 = 1.5 - dist 230 continue if (i1 .EQ. 3) then i1p = i1 + im else i1p = i1 endif if (incc .EQ. 1) go to 260 !----------------------------------------------------------------------- ! j=2 for fine points to the north !----------------------------------------------------------------------- do i = i1p,i2,incc iii = i + ii - 1 psci(iii) = pspc (i-im)*d2 + pspc (i+ip)*d1 sqkci(iii) = sqk0c(i-im)*d2 + sqk0c(i+ip)*d1 do k = 1,kmax uci(k,iii) = upc(k,i-im)*d2 + upc(k,i+ip)*d1 vci(k,iii) = vpc(k,i-im)*d2 + vpc(k,i+ip)*d1 tci(k,iii) = tpc(k,i-im)*d2 + tpc(k,i+ip)*d1 rci(k,iii) = rpc(k,i-im)*d2 + rpc(k,i+ip)*d1 enddo enddo if (im .EQ. 0) go to 290 if (i1 .EQ. 3) then i = i1 iii = i + ii - 1 im = 1 psci (iii) = pspc (i-im)*d2 + pspc (i+ip)*d1 sqkci(iii) = sqk0c(i-im)*d2 + sqk0c(i+ip)*d1 do k = 1,kmax uci(k,iii) = upc(k,i-im)*d2 + upc(k,i+ip)*d1 vci(k,iii) = vpc(k,i-im)*d2 + vpc(k,i+ip)*d1 tci(k,iii) = tpc(k,i-im)*d2 + tpc(k,i+ip)*d1 rci(k,iii) = rpc(k,i-im)*d2 + rpc(k,i+ip)*d1 enddo endif go to 290 !----------------------------------------------------------------------- ! j=jmax-2 for coarse points to the north !----------------------------------------------------------------------- 260 continue do i = i1p,i2,incn iii = i + ii - 1 psni (iii) = pspn (i-im)*d2 + pspn (i+ip)*d1 sqkni(iii) = sqk0n(i-im)*d2 + sqk0n(i+ip)*d1 do k = 1,kmax uni(k,iii) = upn(k,i-im)*d2 + upn(k,i+ip)*d1 vni(k,iii) = vpn(k,i-im)*d2 + vpn(k,i+ip)*d1 tni(k,iii) = tpn(k,i-im)*d2 + tpn(k,i+ip)*d1 rni(k,iii) = rpn(k,i-im)*d2 + rpn(k,i+ip)*d1 enddo enddo if (im .EQ. 0) go to 290 if (i1 .EQ. 3) then i = i1 iii = i + ii - 1 im = 1 psni (iii) = pspn (i-im)*d2 + pspn (i+ip)*d1 sqkni(iii) = sqk0n(i-im)*d2 + sqk0n(i+ip)*d1 do k = 1,kmax uni(k,iii) = upn(k,i-im)*d2 + upn(k,i+ip)*d1 vni(k,iii) = vpn(k,i-im)*d2 + vpn(k,i+ip)*d1 tni(k,iii) = tpn(k,i-im)*d2 + tpn(k,i+ip)*d1 rni(k,iii) = rpn(k,i-im)*d2 + rpn(k,i+ip)*d1 enddo endif 290 continue enddo endif !----------------------------------------------------------------------- ! fluxes for inner nest for j=3 or j=jmax-1 !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! added to decouple the south and north fluxes !----------------------------------------------------------------------- if (incc .NE. incs) then do ii = 1,insc1 dist = FLOAT(2*ii - 1)/FLOAT(2*insc1) if (ii .GT. insc2) go to 300 ip = 0 im = insc1 d1 = .5 + dist d2 = .5 - dist go to 310 300 continue ip = insc1 im = 0 d1 = -.5 + dist d2 = 1.5 - dist 310 continue if (i1 .EQ. 3) then i1p = i1 + im else i1p = i1 endif if (incc .EQ. 1) go to 340 !----------------------------------------------------------------------- ! j=jmax-1 for fine points to the south !----------------------------------------------------------------------- do i = i1p,i2,incc iii = i + ii - 1 psci (iii) = pspc (i-im)*d2 + pspc (i+ip)*d1 sqkci(iii) = sqk0c(i-im)*d2 + sqk0c(i+ip)*d1 do k = 1,kmax uci(k,iii) = upc(k,i-im)*d2 + upc(k,i+ip)*d1 vci(k,iii) = vpc(k,i-im)*d2 + vpc(k,i+ip)*d1 tci(k,iii) = tpc(k,i-im)*d2 + tpc(k,i+ip)*d1 rci(k,iii) = rpc(k,i-im)*d2 + rpc(k,i+ip)*d1 enddo enddo if (im .EQ. 0) go to 370 if (i1 .EQ. 3) then i = i1 iii = i + ii - 1 im = 1 psci (iii) = pspc (i-im)*d2 + pspc (i+ip)*d1 sqkci(iii) = sqk0c(i-im)*d2 + sqk0c(i+ip)*d1 do k = 1,kmax uci(k,iii) = upc(k,i-im)*d2 + upc(k,i+ip)*d1 vci(k,iii) = vpc(k,i-im)*d2 + vpc(k,i+ip)*d1 tci(k,iii) = tpc(k,i-im)*d2 + tpc(k,i+ip)*d1 rci(k,iii) = rpc(k,i-im)*d2 + rpc(k,i+ip)*d1 enddo endif go to 370 !----------------------------------------------------------------------- ! j=3 for coarse points to the south !----------------------------------------------------------------------- 340 continue do i = i1p,i2,incs iii = i + ii - 1 pssi (iii) = psps (i-im)*d2 + psps (i+ip)*d1 sqksi(iii) = sqk0s(i-im)*d2 + sqk0s(i+ip)*d1 do k = 1,kmax usi(k,iii) = ups(k,i-im)*d2 + ups(k,i+ip)*d1 vsi(k,iii) = vps(k,i-im)*d2 + vps(k,i+ip)*d1 tsi(k,iii) = tps(k,i-im)*d2 + tps(k,i+ip)*d1 rsi(k,iii) = rps(k,i-im)*d2 + rps(k,i+ip)*d1 enddo enddo if (im .EQ. 0) go to 370 if (i1 .EQ. 3) then i = i1 iii = i + ii - 1 im = 1 pssi (iii) = psps (i-im)*d2 + psps (i+ip)*d1 sqksi(iii) = sqk0s(i-im)*d2 + sqk0s(i+ip)*d1 do k = 1,kmax usi(k,iii) = ups(k,i-im)*d2 + ups(k,i+ip)*d1 vsi(k,iii) = vps(k,i-im)*d2 + vps(k,i+ip)*d1 tsi(k,iii) = tps(k,i-im)*d2 + tps(k,i+ip)*d1 rsi(k,iii) = rps(k,i-im)*d2 + rps(k,i+ip)*d1 enddo endif 370 continue enddo endif !----------------------------------------------------------------------- ! in fine nest this will take care of all the rows except for the ! j=2 and j=jmax-1 rows !----------------------------------------------------------------------- if (incc .EQ. incr .AND. incc .EQ. insc) then do i = i1,i2 psci (i) = pspc (i) sqkci(i) = sqk0c(i) enddo do i = i1,i2 do k = 1,kmax uci(k,i) = upc(k,i) vci(k,i) = vpc(k,i) tci(k,i) = tpc(k,i) rci(k,i) = rpc(k,i) enddo enddo endif !----------------------------------------------------------------------- ! for the j=2 and j=jmax-1 rows in the fine nest, we must define an ! array with the original coarse data to obtain the flux on the ! coarse to coarse side of the row. for all other rows, these arrays ! will be the original "xci" arrays !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! j=2 row of the fine mest !----------------------------------------------------------------------- if (j .EQ. 2 .AND. incc .NE. 1) then do i = i1,i2 pscsi (i) = pspc (i) sqkcsi(i) = sqk0c(i) enddo do i = i1,i2 do k = 1,kmax ucsi(k,i) = upc(k,i) vcsi(k,i) = vpc(k,i) tcsi(k,i) = tpc(k,i) rcsi(k,i) = rpc(k,i) enddo enddo else do i = i1,i2 pscsi (i) = psci (i) sqkcsi(i) = sqkci(i) enddo do i = i1,i2 do k = 1,kmax ucsi(k,i) = uci(k,i) vcsi(k,i) = vci(k,i) tcsi(k,i) = tci(k,i) rcsi(k,i) = rci(k,i) enddo enddo endif !----------------------------------------------------------------------- ! j=jmax-1 row of the fine mesh !----------------------------------------------------------------------- if (j .EQ. jmax - 1 .AND. incc .NE. 1) then do i = i1,i2 pscni (i) = pspc(i) enddo do i = i1,i2 do k = 1,kmax ucni(k,i) = upc(k,i) vcni(k,i) = vpc(k,i) tcni(k,i) = tpc(k,i) rcni(k,i) = rpc(k,i) enddo enddo else do i = i1,i2 pscni (i) = psci (i) enddo do i = i1,i2 do k = 1,kmax ucni(k,i) = uci(k,i) vcni(k,i) = vci(k,i) tcni(k,i) = tci(k,i) rcni(k,i) = rci(k,i) enddo enddo endif !----------------------------------------------------------------------- ! all rows except the j=jmax-2 row !----------------------------------------------------------------------- if (incn .EQ. incr) then do i = i1,i2 psni (i) = pspn (i) sqkni(i) = sqk0n(i) enddo do i = i1,i2 do k = 1,kmax uni(k,i) = upn(k,i) vni(k,i) = vpn(k,i) tni(k,i) = tpn(k,i) rni(k,i) = rpn(k,i) enddo enddo endif !----------------------------------------------------------------------- ! all rows except the j=3 row !----------------------------------------------------------------------- if (incs .EQ. insc) then do i = i1,i2 pssi (i) = psps (i) sqksi(i) = sqk0s(i) enddo do i = i1,i2 do k = 1,kmax usi(k,i) = ups(k,i) vsi(k,i) = vps(k,i) tsi(k,i) = tps(k,i) rsi(k,i) = rps(k,i) enddo enddo endif !----------------------------------------------------------------------- ! compute averaged southern fluxes interface values at northern ! interface of previous row are saved as southern interface values ! of current row !----------------------------------------------------------------------- if (j .EQ. 1) go to 590 do i = i1,i2 do k = 1,kmax dus1(k,i) = rlyss(k,i)*csrlts(k,i)*(ucsi(k,i)/csrlt - & usi(k,i)/csrlss) dvs1(k,i) = rlyss(k,i)*(vcsi(k,i) - vsi(k,i)) dts1(k,i) = rlyss(k,i)*(tcsi(k,i) - tsi(k,i)) drs1(k,i) = rlyss(k,i)*(rcsi(k,i) - rsi(k,i)) enddo enddo do i = i1,i2 do k = 1,kmax dfs1(k,i) = SQRT(4.*dvs1(k,i)**2 + dus1(k,i)**2) enddo enddo do i = i1,i2 ely = 0.5*(dlatpc(i) + dlatps(i))*arad psb = (pssi (i)*(1. - gamas(i)) + pscsi (i)*gamas(i)) sqkb = (sqksi(i)*(1. - gamas(i)) + sqkcsi(i)*gamas(i)) coys1(i) = psb*sqkb*ely*ely enddo do i = i1,i2 do k = 1,kmax reys1(k,i) = coys1(i) enddo enddo do i = i1,i2 do k = 1,kmax sxys1(k,i) = reys1(k,i)*dfs1(k,i)*dus1(k,i) syys1(k,i) = reys1(k,i)*dfs1(k,i)*dvs1(k,i)*2 htys1(k,i) = reys1(k,i)*dfs1(k,i)*dts1(k,i) hrys1(k,i) = reys1(k,i)*dfs1(k,i)*drs1(k,i) enddo enddo dincc = 1./FLOAT(incc) if (incc .EQ. incs) dincc = 1. do i = i1,i2 do k = 1,kmax sxys(k,i) = sxys1(k,i)*dincc syys(k,i) = syys1(k,i)*dincc htys(k,i) = htys1(k,i)*dincc hrys(k,i) = hrys1(k,i)*dincc enddo enddo if (incc .LE. incs) go to 590 do ii = 2,incc ih = ii - 1 do i = i1,i2 do k = 1,kmax sxys(k,i) = sxys(k,i) + sxys1(k,i+ih)*dincc syys(k,i) = syys(k,i) + syys1(k,i+ih)*dincc htys(k,i) = htys(k,i) + htys1(k,i+ih)*dincc hrys(k,i) = hrys(k,i) + hrys1(k,i+ih)*dincc enddo enddo enddo !----------------------------------------------------------------------- ! calculate fluxes at northern interface !----------------------------------------------------------------------- 590 continue if (j .EQ. jmax .AND. nstflg .EQ. 1) go to 710 if (j .EQ. jmax .AND. nstflg .EQ. 0) go to 670 do i = i1,i2 do k = 1,kmax dun1(k,i) = rlynn(k,i)*csrltn(k,i)*(uni(k,i)/csrlnn - & ucni(k,i)/csrlt) dvn1(k,i) = rlynn(k,i)*(vni(k,i) - vcni(k,i)) dtn1(k,i) = rlynn(k,i)*(tni(k,i) - tcni(k,i)) drn1(k,i) = rlynn(k,i)*(rni(k,i) - rcni(k,i)) enddo enddo do i = i1,i2 do k = 1,kmax dfn1(k,i) = SQRT(4.*dvn1(k,i)**2 + dun1(k,i)**2) enddo enddo do i = i1,i2 ely = 0.5*(dlatpc(i) + dlatpn(i))*arad psb = (psni (i)*(1. - gaman(i)) + pscni(i)*gaman(i)) sqkb = (sqkni(i)*(1. - gaman(i)) + sqkci(i)*gaman(i)) coyn1(i) = psb*sqkb*ely*ely enddo do i = i1,i2 do k = 1,kmax reyn1(k,i) = coyn1(i) enddo enddo do i = i1,i2 do k = 1,kmax sxyn1(k,i) = reyn1(k,i)*dfn1(k,i)*dun1(k,i) syyn1(k,i) = reyn1(k,i)*dfn1(k,i)*dvn1(k,i)*2 htyn1(k,i) = reyn1(k,i)*dfn1(k,i)*dtn1(k,i) hryn1(k,i) = reyn1(k,i)*dfn1(k,i)*drn1(k,i) enddo enddo !----------------------------------------------------------------------- ! average fluxes for northern interface !----------------------------------------------------------------------- dincc = 1./FLOAT(incc) if (incc .EQ. incn) dincc = 1. do i = i1,i2 do k = 1,kmax sxyn(k,i) = sxyn1(k,i)*dincc syyn(k,i) = syyn1(k,i)*dincc htyn(k,i) = htyn1(k,i)*dincc hryn(k,i) = hryn1(k,i)*dincc enddo enddo if (incc .LE. incn) go to 670 do ii = 2,incc ih = ii - 1 do i = i1,i2 do k = 1,kmax sxyn(k,i) = sxyn(k,i) + sxyn1(k,i+ih)*dincc syyn(k,i) = syyn(k,i) + syyn1(k,i+ih)*dincc htyn(k,i) = htyn(k,i) + htyn1(k,i+ih)*dincc hryn(k,i) = hryn(k,i) + hryn1(k,i+ih)*dincc enddo enddo enddo !----------------------------------------------------------------------- ! boundary conditions on northern and southern fluxes !----------------------------------------------------------------------- 670 continue if (j .EQ. jmax) go to 730 if (j .GT. 1) go to 750 if (nstflg .EQ. 1) go to 690 !----------------------------------------------------------------------- ! bottom southernmost row of coarsest grid !----------------------------------------------------------------------- do i = i1,i2 do k = 1,kmax sxys(k,i) = sxyn(k,i) syys(k,i) = syyn(k,i) htys(k,i) = htyn(k,i) hrys(k,i) = hryn(k,i) enddo enddo go to 750 !----------------------------------------------------------------------- ! bottom southernmost row of current nested grid !----------------------------------------------------------------------- 690 continue ii1 = i1 if (i1 .EQ. imax-1) ii1 = (imax - 4)/ngr + 3 ii2 = (i2 - i1 + 1)/incc + ii1 - 1 i = i1 - incc iii = 0 do iir = ii1,ii2 iii = iii + 1 ii = (ii1 - 3 + incc)/incc + iii + 1 i = i + incc do k = 1,kmax !----------------------------------------------------------------------- ! sxys(k,i) = bnfxm2(k,1,ii,17, ngd-1) ! syys(k,i) = bnfxm2(k,1,ii,18, ngd-1) ! htys(k,i) = bnfxm2(k,1,ii,19,ngd-1) ! hrys(k,i) = bnfxm2(k,1,ii,20,ngd-1) !----------------------------------------------------------------------- sxys(k,i) = bnfxm2(k,1,ii,23, ngd-1) syys(k,i) = bnfxm2(k,1,ii,24, ngd-1) htys(k,i) = bnfxm2(k,1,ii,25,ngd-1) hrys(k,i) = bnfxm2(k,1,ii,26,ngd-1) enddo enddo go to 750 !----------------------------------------------------------------------- ! top northernmost row of current nested grid !----------------------------------------------------------------------- 710 continue ii1 = i1 if (i1 .EQ. imax-1) ii1 = (imax - 4)/ngr + 3 ii2 = (i2 - i1 + 1)/incc + ii1 - 1 i = i1 - incc iii = 0 do iir = ii1,ii2 iii = iii + 1 ii = (ii1 - 3 + incc)/incc + iii + 1 i = i + incc do k = 1,kmax !----------------------------------------------------------------------- ! sxyn(k,i) = bnfxm2(k,2,ii,17 ,ngd-1) ! syyn(k,i) = bnfxm2(k,2,ii,18 ,ngd-1) ! htyn(k,i) = bnfxm2(k,2,ii,19,ngd-1) ! hryn(k,i) = bnfxm2(k,2,ii,20,ngd-1) !----------------------------------------------------------------------- sxyn(k,i) = bnfxm2(k,2,ii,23 ,ngd-1) syyn(k,i) = bnfxm2(k,2,ii,24 ,ngd-1) htyn(k,i) = bnfxm2(k,2,ii,25,ngd-1) hryn(k,i) = bnfxm2(k,2,ii,26,ngd-1) enddo enddo go to 750 !----------------------------------------------------------------------- ! top row of coarsest grid computed above with average fluxes for ! northern interface. outer boundary flux (at jmax) = flux at ! jmax-1 !----------------------------------------------------------------------- 730 continue !----------------------------------------------------------------------- ! top northernmost row of current outside grid !----------------------------------------------------------------------- do i = i1,i2 do k = 1,kmax sxyn(k,i) = sxys(k,i) syyn(k,i) = syys(k,i) htyn(k,i) = htys(k,i) hryn(k,i) = hrys(k,i) enddo enddo 750 continue i1p = i1 + 1 i2p = i2 + 1 i11 = i1 i22 = i2 + 1 if (i1 .EQ. imax) then i22 = i2 endif if (i1 .EQ. 1) then i11 = i1 + 1 endif if (nstflg .EQ. 1 .AND. i1 .EQ. 3) then i11 = i1 + 1 endif do i = i11,i22 gamaw(i) = 0.5 rlxw(i) = 1./(.5*(dlonpc(i) + dlonpc(i-incc))*arad* & COS(rlatpc(i))) enddo !----------------------------------------------------------------------- ! compute western fluxes !----------------------------------------------------------------------- i11 = i1 i22 = i2 + 1 if(i1 .EQ. 1. AND. nstflg .EQ. 0) then i11 = i1 + 1 i22 = i2 + 1 endif if(i1 .EQ. 1. AND. nstflg .EQ. 1) then i11 = i1 + 1 i22 = i2 endif if(i2 .EQ. imax. AND. nstflg .EQ. 0) then i22 = i2 i11 = i1 endif if(i2 .EQ. imax. AND. nstflg .EQ. 1) then i11 = i1 + 1 endif if ((i1 .EQ. 3 .OR. i1 .EQ. 1) .AND. nstflg .EQ. 1) then i11 = i1 + 1 endif if ((i2 .EQ. imax-2 .OR. i2 .EQ. imax) .AND. nstflg .EQ. 1) then i22 = i2 endif do i = i11,i22 psbarw(i) = pspc (i)*(1. - gamaw(i)) + pspc (i-incc)*gamaw(i) sqkw (i) = sqk0c(i)*(1. - gamaw(i)) + sqk0c(i-incc)*gamaw(i) enddo do i = i11,i22 do k = 1,kmax duw(k,i) = rlxw(i)*(upc(k,i) - upc(k,i-incc)) - tnrlt* & (vpc(k,i-incc)*gamaw(i) + vpc(k,i)* & (1. - gamaw(i)))/arad dvw(k,i) = rlxw(i)*(vpc(k,i) - vpc(k,i-incc)) dtw(k,i) = rlxw(i)*(tpc(k,i) - tpc(k,i-incc)) drw(k,i) = rlxw(i)*(rpc(k,i) - rpc(k,i-incc)) enddo enddo do i = i11,i22 coxw(i) = 0.5*(dlonpc(i) + dlonpc(i-incc))*arad*COS(rlatpc(i)) enddo if (nstflg .EQ. 1 .AND. i1 .EQ. 3) then coxw(i1) = 0.5*(dlonpc(i1) + dlonpc(i1 - 1))* & arad*COS(rlatpc(i1)) endif if (nstflg .EQ. 1 .AND. i2 .EQ. imax - 2) then coxw(i2p) = 0.5*(dlonpc(i2p) + dlonpc(i2p - incc))*arad* & COS(rlatpc(i2p - incc)) endif if (nstflg .EQ. 1 .AND. i1 .EQ. 3) then i11 = 3 endif if (nstflg .EQ. 1 .AND. i2 .EQ. imax-2) then i22 = i2 + 1 endif do i = i11,i22 do k = 1,kmax dfw(k,i) = SQRT(4.*duw(k,i)**2 + dvw(k,i)**2) enddo enddo do i = i11,i22 do k = 1,kmax dfw(k,i) = dfw(k,i)*coxw(i)*coxw(i)*psbarw(i)*sqkw(i) enddo enddo do i = i11,i22 do k = 1,kmax sxxw(k,i) = dfw(k,i)*duw(k,i)*2. syxw(k,i) = dfw(k,i)*dvw(k,i) htxw(k,i) = dfw(k,i)*dtw(k,i) hrxw(k,i) = dfw(k,i)*drw(k,i) enddo enddo !----------------------------------------------------------------------- ! outermost grid east and west boundary fluxes at i=1 and imax+1 ! are set to fluxes at i=2 and imax !----------------------------------------------------------------------- if (i1 .EQ. 1 .AND. nstflg .EQ. 0) then do k = 1,kmax sxxw(k,1) = sxxw(k,2) syxw(k,1) = syxw(k,2) htxw(k,1) = htxw(k,2) hrxw(k,1) = hrxw(k,2) enddo endif if (i2 .EQ. imax .AND. nstflg .EQ. 0) then do k = 1,kmax sxxw(k,imax + 1) = sxxw(k,imax) syxw(k,imax + 1) = syxw(k,imax) htxw(k,imax + 1) = htxw(k,imax) hrxw(k,imax + 1) = hrxw(k,imax) enddo endif !----------------------------------------------------------------------- ! store boundary fluxes at edge of inner nest hole for dynamical ! interface of next finer grid !----------------------------------------------------------------------- if (ngd .NE. nnst) then if (j .LT. js .OR. j .GE. jn .OR. i2 .LE. iw .OR. & i1 .GT. ie) go to 880 if (j .NE. js) go to 830 jf = 1 go to 840 830 continue if (j .NE. jn - 1) go to 860 jf = 2 840 continue if (i1 .EQ. ie) go to 860 if (i1 .GE. iw + 1 .AND. i1 .LE. ie - 1) then iwp = i1 else iwp = iw + 1 endif if (i2 .GE. iw + 1 .AND. i2 .LE. ie - 1) then iem = i2 else iem = ie - 1 endif do i = iwp,iem ii = i - iw do k = 1,kmax !----------------------------------------------------------------------- ! bnfxm2(k,jf,ii,17, ngd) = sxyn(k,i) ! bnfxm2(k,jf,ii,18, ngd) = syyn(k,i) ! bnfxm2(k,jf,ii,19,ngd) = htyn(k,i) ! bnfxm2(k,jf,ii,20,ngd) = hryn(k,i) !----------------------------------------------------------------------- bnfxm2(k,jf,ii,23, ngd) = sxyn(k,i) bnfxm2(k,jf,ii,24, ngd) = syyn(k,i) bnfxm2(k,jf,ii,25,ngd) = htyn(k,i) bnfxm2(k,jf,ii,26,ngd) = hryn(k,i) enddo enddo 860 continue if (j .EQ. js) go to 880 jf = j - js do i = i1,i2 if (i .EQ. iw + 1) then ii = 1 endif if (i .EQ. ie) then ii = 2 endif if (i .EQ. iw + 1 .OR. i .EQ. ie) then do k = 1,kmax !----------------------------------------------------------------------- ! bnfxm2(k,ii,jf,21, ngd) = sxxw(k,i) ! bnfxm2(k,ii,jf,22, ngd) = syxw(k,i) ! bnfxm2(k,ii,jf,23,ngd) = htxw(k,i) ! bnfxm2(k,ii,jf,24,ngd) = hrxw(k,i) !----------------------------------------------------------------------- bnfxm2(k,ii,jf,27, ngd) = sxxw(k,i) bnfxm2(k,ii,jf,28, ngd) = syxw(k,i) bnfxm2(k,ii,jf,29,ngd) = htxw(k,i) bnfxm2(k,ii,jf,30,ngd) = hrxw(k,i) enddo endif enddo endif !----------------------------------------------------------------------- ! compute diffusive flux divergence !----------------------------------------------------------------------- 880 continue if (nstflg .EQ. 1) go to 900 if (j .NE. 1 .AND. j .NE. jmax) go to 900 do i = i1,i2 do k = 1,kmax hfu(k,i) = (sxxw(k,i+incc) - sxxw(k,i))*wtw(k,i) hfv(k,i) = (syxw(k,i+incc) - syxw(k,i))*wtw(k,i) hft(k,i) = (htxw(k,i+incc) - htxw(k,i))*wtw(k,i) hfr(k,i) = (hrxw(k,i+incc) - hrxw(k,i))*wtw(k,i) enddo enddo go to 920 900 continue do i = i1,i2 do k = 1,kmax hfu(k,i) = (sxxw(k,i+incc) - sxxw(k,i))*wtw(k,i) + sxyn(k,i)* & wtn(k,i) - sxys(k,i)*wts(k,i) hfv(k,i) = (syxw(k,i+incc) - syxw(k,i))*wtw(k,i) + syyn(k,i)* & wtn(k,i) - syys(k,i)*wts(k,i) hft(k,i) = (htxw(k,i+incc) - htxw(k,i))*wtw(k,i) + htyn(k,i)* & wtn(k,i) - htys(k,i)*wts(k,i) hfr(k,i) = (hrxw(k,i+incc) - hrxw(k,i))*wtw(k,i) + hryn(k,i)* & wtn(k,i) - hrys(k,i)*wts(k,i) enddo enddo 920 continue if (nstflg .EQ. 1) then if (i1 .EQ. 3) then dyinc = dlatpc(i1)/dlatpc(i1-1) cumdy(1,jj) = cumdy(1,jj) + dyinc do k = 1,kmax cumhid(k,1,jj,1) = cumhid(k,1,jj,1) + sxxw(k,i1)*dyinc cumhid(k,1,jj,2) = cumhid(k,1,jj,2) + syxw(k,i1)*dyinc cumhid(k,1,jj,3) = cumhid(k,1,jj,3) + htxw(k,i1)*dyinc cumhid(k,1,jj,4) = cumhid(k,1,jj,4) + hrxw(k,i1)*dyinc enddo if (jrow .EQ. njrow) then if (ABS(cumdy(1,jj) - 1.) .GT. .1) then write(stdout,930) ngd, j, (cumdy(i,jj), i=1,2) write(stdout,935) dlatpc(i1), dlatpc(i1-1), i1, & dlatpc(i1) 930 FORMAT(1X, ' GRID = ', I5, ' ROW = ', I5, ' BAD CUMDY I1 = ', & 2E20.5) 935 FORMAT(1X, ' DEL = ', 2E20.5, 2X, I5, E20.5) stop endif endif endif if (i2 .EQ. imax - 2) then dyinc = dlatpc(i2-ngr+1)/dlatpc(i2+1) cumdy(2,jj) = cumdy(2,jj) + dyinc do k = 1,kmax cumhid(k,2,jj,1) = cumhid(k,2,jj,1) + sxxw(k,i2p)*dyinc cumhid(k,2,jj,2) = cumhid(k,2,jj,2) + syxw(k,i2p)*dyinc cumhid(k,2,jj,3) = cumhid(k,2,jj,3) + htxw(k,i2p)*dyinc cumhid(k,2,jj,4) = cumhid(k,2,jj,4) + hrxw(k,i2p)*dyinc enddo if (jrow .EQ. njrow) then if (ABS(cumdy(2,jj) - 1.) .GT. .1) then write(stdout,940) ngd, j, (cumdy(i,jj), i=1,2) write(stdout,945) dlatpc(i2-ngr+1),dlatpc(i2+1),i2, & dlatpc(i2) 940 FORMAT(1X, ' GRID = ', I5, ' ROW = ', I5, ' BAD CUMDY I2 = ', & 2E20.5) 945 FORMAT(1X, ' DELX = ', 2E20.5, 2X, I5, E20.5) stop endif endif endif endif return end subroutine HDIFSN