c Extracted the COMMON block /GIANT4/ to a header file. The array PHNF c is transferred to a new COMMON block /GIANT5/. c SJP 2009/03/12 c c Parallel compiler directives replaced with equivalent OpenMP instructions. c SJP 2001/12/12 c c $Log: dynmnl.f,v $ c Revision 1.28 2001/06/04 02:27:00 rot032 c Changes from LDR to bring sulfur cycle to run P41 c c Revision 1.27 2001/02/28 04:36:40 rot032 c Further tidy ups from HBG c c Revision 1.26 2001/02/22 05:34:43 rot032 c Changes from HBG to complete concatenation of NH/SH latitudes c c Revision 1.25 1999/06/22 04:11:12 rot032 c Add geopotential height diagnostic and tidy datard.f. c c Revision 1.24 1999/06/16 06:21:54 rot032 c HBG changes to V5-3 c c Revision 1.23 1999/05/20 06:23:55 rot032 c HBG changes to V5-2 c c Revision 1.22 1997/12/17 23:22:55 ldr c Changes from MRD for parallel execution on NEC. c c Revision 1.21 1996/10/24 01:02:39 ldr c Comments, Legendre transforms and tidy-ups from TIE c c Revision 1.20 1996/03/21 03:18:36 ldr c Changes from TIE to remove nsemilag.eq.0 and tidy physical constants c c Revision 1.19 1995/08/31 04:30:41 ldr c Tidy-ups from HBG affecting character arrays and generality of NL c c Revision 1.18 1995/08/18 05:52:51 ldr c HBG's changes to V4-7-12l for hybrid coordinate. c c Revision 1.17 1995/07/05 06:14:32 ldr c Changes from Dave Micklethwaite for Cray parallel execution tidied by LDR c and merged into V4-7-2l. c c Revision 1.15.1.1 1995/07/03 05:48:41 ldr c Mickles' changes to V4-5-30l for Cray multiprocessing, tidied by LDR. c c Revision 1.16 1994/08/08 17:21:05 ldr c Strip off excessive RCS comments at top of file. c c Revision 1.15 93/12/17 15:32:12 ldr c Hack V4-4-52l to change all continuation chars to & c c Revision 1.14 93/10/05 13:05:57 ldr c Changes to V4-4-15l from HBG for T63 and coupled model c c Revision 1.13 93/08/03 11:25:38 ldr c Replace intrinsic functions with generic forms where possible. c c Revision 1.12 93/07/14 14:54:13 ldr c The linear terms generated by P*f are now included in the nonlinear terms c (ready for implicit vorticity)(HBG). c c Revision 1.11 93/06/16 14:11:10 ldr c HBG changes for Chen temperature variable. c c Revision 1.10 93/03/15 15:45:11 ldr c HBG's changes to combine uvpgd and sltjmcg common blocks. c c Revision 1.9 93/01/26 17:34:17 ldr c Removed factor 100 from evaluation of dpsdt. c c Revision 1.8 92/12/10 09:56:18 ldr c Replace np's with nlp's for compatibility with ogcm. c c INPUT/OUTPUT: c Input: from common/cnsta in CNSTA.f c dsk - sigma thicknesses (1 to nl) c sig - sigma values sigk - sigma**cappa c c from common/const in CONST.f c cappa - specific gas const dry air/spec heat dry air at c constant pressure c tomg - 2 X omega c c from common/fewflags in FEWFLAGS.f c chenflag - flag for CHEN version of model c hybrid - if true, use new hybrid vertical coordinate c rather than sigma coordinate c c from common/gausl in GAUSL.f c acsq - 1.0/(sia(lg)*sia(lg)) c coa - sin latitude c sia - cosine latitude c c from common/hybarr in HYBARR.f c anf, bnf, bnh, dadnf, dbdnf - hybrid vertical coordinate c arrays, f=full level, h=half levels c c from common/msp21 in this subroutine c psrk - chen surface pressure variable c c from common/si in SI.f c am - geopotential height matrix c amh- " " " , half levels c c from common/worka in WORKA.f c emean - preset global mean energy const/model level c phisb - global mean geopotential height c tmean - preset isothermal temperature c c from common/work1x in this subroutine c pn - latitude row of surface pressure c cpn - gradients of pn c rmn - surface pressure weighted moisture c ten - surface pressure weighted grid point temp c un - zonal wind velocity m/s c vn - meridional wind velocity m/s c von - grid point vorticity c dvn - grid point divergence c pln - D(P*)/D(LONG c c Output: from common/giant1 in this subroutine c aon, bon, eon, fon, gon, hon - grid values of A, B, E, c F, G, H before Fast Fourier Transforms c c from common/prtar in this subroutine c dvgg - global array of divergence c c from common/uvpgd in UVPGD.f c dpsdt - omega/sigma at full levels c sdot - d(sigma)/d(t) c c In/Out: from common/dynav in this subroutine c ps - surface pressure ww - omega c t - zonally averaged temp at all levels c u - zonal wind v - meridional wind, ww c c from common/giant4 in this subroutine c psl - sea level pressure c rpn - pure moisture c tdn - non surface pressure weighted temp. c tpn - pure temp(real temp) c upn - u*cos(phi)/erad vpn - v*cos(phi)/erad c z4 - grid elevations [non-spectral] c ron - measure of east-west stress at surface c son - measure of north-south stress at surface c c from common/hybrpr in HYBRPR.f c muf - pressure derivative with respect to coordinates, at c full levels c prf - pressure at full levels c c from common/vertv in VERTV.f c pssd, psv, pvbqb, pvbtb, pvbub, pvqb, pvtb, pvub c - zonally averaged meridional fluxes c subroutine dynmnl(lg,lstat,zonpsl,lgn) implicit none !$OMP THREADPRIVATE ( /GIANT1/ ) !$OMP THREADPRIVATE ( /GIANT4/ ) !$OMP THREADPRIVATE ( /GIANT5/ ) !$OMP THREADPRIVATE ( /HYBRPR/ ) !$OMP THREADPRIVATE ( /WORK1X/ ) C Global parameters include 'PARAMS.f' include 'PHYSPARAMS.f' real qcrit parameter (qcrit=2.e-6) C Argument list integer lg logical lstat real zonpsl(2) integer lgn C Local shared common blocks (see TASK COMMON/TASKLOCAL above) real aon,bon,eon,fon,gon,hon common/giant1/aon(ln2,nl) & ,bon(ln2,nl),eon(ln2,nl),fon(ln2,nl),gon(ln2,nl),hon(ln2,nl) include 'GIANT4.f' include 'GIANT5.f' include 'HYBRPR.f' real un,vn,pn,dvn,von,ten,cpn,pln,rmn common/work1x/un(ln2,nl),vn(ln2,nl),pn(ln2),dvn(ln2,nl) & ,von(ln2,nl),ten(ln2,nl),cpn(ln2),pln(ln2),rmn(ln2,nl) C Global data blocks include 'CNSTA.f' include 'FEWFLAGS.f' include 'GAUSL.f' include 'HYBARR.f' include 'SI.f' include 'UVPGD.f' include 'VERTV.f' !Coupled model include 'WORKA.f' real u,v,t,ww,ps common/dynav/u(lat,2,nl),v(lat,2,nl),t(lat,2,nl) &,ww(lat,2,nlm),ps(lat,2) real psrk common/msp21/psrk(ln2,lat) real devap,drain,dvgg common/prtar/devap(ln2,lat),drain(ln2,lat),dvgg(ln2,lat) c real omegac(ln2,nl,lat) c common/omblk/omegac C Local work arrays and variables real dvhn(ln2,nlp),wbn(ln2,nlp),uh(ln2,nlp),vh(ln2,nlp) & ,th(ln2,nlp) real phist(ln2),phn(ln2,nl) & ,dvfn(ln2,nl),phnh(ln2,nlp),tpnv(ln2,nl),phnv(ln2,nl) real temu(ln2,nl),temv(ln2,nl),temt(ln2,nl) real pnk(ln2) real fff1(ln2,nl),fff2(ln2,nl),fff3(ln2,nl),phnhat(ln2,nl) real txn1(ln2,nl),txn2(ln2,nl),txn3(ln2,nl) real plk(ln2,2),tlk(ln2,2) real dpsdtk(ln2,nl) real tomgns(ln2) integer j integer k integer ma integer mg integer ns real c real con real conr real csqlg real delt2 real dskt real etem real fchen real fhyb real flong real psvz real psz real pvqbz real pvtbz real pvubz real pxlev real p00 real qz real siai real tem real temz1 real txlev real tz real t0ec real t00ec real t1ec real upnk real uz real vonk real vpnk real vz real wsx real wzk real xx1 real xx2 C Local data, functions etc C Start code : ---------------------------------------------------------- CHEN t00ec=288.0 t1ec=0.6652*t00ec t0ec=t00ec-t1ec p00=1000.0 CHEN csqlg=acsq(lg) siai=erad/sia(lg) C**** COMPUTE (muf=dp/dn) WEIGHTS, PRESSURE (full levels) do 80 k=1,nl do 80 mg=1,ln2 muf(mg,k)=dadnf(k)+(dbdnf(k)*pn(mg)) 80 prf(mg,k)=anf(k)+(bnf(k)*pn(mg)) C* MG SIGNIFIES LONGITUDE GRID POINT W TO E C**** C**** NON-LINEAR DYNAMICS C**** C* COMPUTE VERT INTEGRAL OF DIVERGENCE TO ALL HALF LEVS DVHN() c* and to all full levels dvfn() do 201 mg=1,ln2 201 dvhn(mg,nlp)=0.0 do 120 j=1,nl k=nlp-j do 120 mg=1,ln2 dvhn(mg,k)=dvhn(mg,k+1)+dsk(k)*dvn(mg,k) 120 dvfn(mg,k)=dvhn(mg,k+1)+0.5*dsk(k)*dvn(mg,k) C* COMPUTE VERTICAL VELOCITIES do 1211 mg=1,ln2 wbn(mg,1)=0.0 1211 wbn(mg,nlp)=0.0 do 121 k=2,nl do 121 mg=1,ln2 121 wbn(mg,k)=bnh(k)*dvhn(mg,1)-dvhn(mg,k) C**** COMPUTE NON [muf(mg,k)] WEIGHTED VALUES AT (full levels) do 123 k=1,nl do 123 mg=1,ln2 upn(mg,k)=un(mg,k)/muf(mg,k) vpn(mg,k)=vn(mg,k)/muf(mg,k) 123 tdn(mg,k)=ten(mg,k)/muf(mg,k) If(chenflag)Then CHEN fchen=1.0 do 1235 mg=1,ln2 1235 pnk(mg)=(pn(mg)/p00)**cappa do 1237 k=1,nl do 1237 mg=1,ln2 if(hybrid)then tpn(mg,k)=tdn(mg,k)+t0ec+t1ec*(prf(mg,k)/p00)**cappa else tpn(mg,k)=tdn(mg,k)+t0ec+t1ec*sigk(k)*pnk(mg) endif 1237 continue CHEN Else fchen=0.0 do 1238 k=1,nl do 1238 mg=1,ln2 1238 tpn(mg,k)=tdn(mg,k)+tmean(k) End If do 124 k=1,nl do 124 mg=1,ln2 c rpn(mg,k)=rmn(mg,k)/muf(mg,k) rpn(mg,k)=rmn(mg,k) c Calculate virtual temperature 124 tpnv(mg,k)=tpn(mg,k)*(epsil+rpn(mg,k))/(epsil*(1.+rpn(mg,k))) C* COMPUTE HALF LEVEL VALUES FOR VERT FLUXES do 125 k=2,nl dskt=dsk(k-1)+dsk(k) xx1=dsk(k)/dskt xx2=dsk(k-1)/dskt do 125 mg=1,ln2 uh(mg,k)=0.5*(upn(mg,k-1)+upn(mg,k))*wbn(mg,k) vh(mg,k)=0.5*(vpn(mg,k-1)+vpn(mg,k))*wbn(mg,k) 125 th(mg,k)=(xx1*tdn(mg,k-1)+xx2*tdn(mg,k))*wbn(mg,k) C**** COMPUTE GEOPOTENTIAL TERMS C**** phist() = GEOPOT.PERT. AT EARTH SURFACE) If(chenflag)Then CHEN do 2072 mg=1,ln2 phist(mg)=cp*(t0ec*log(pnk(mg)/psrk(mg,lg))+ & t1ec*(pnk(mg)-psrk(mg,lg))) 2072 continue CHEN Else do 2074 mg=1,ln2 2074 phist(mg)=z4(mg)*grav-phisb End If C**** phnv() = (full level) GEOPOT. DIFFERENCE DUE TO THE DIFFERENCE C**** BETWEEN VIRTUAL- AND ACTUAL TEMPERATURE C**** phn() = (full level) GEOPOT.PERT. DEPENDING UPON T'=tpn(),tmean(). C**** phnh() = (half level) GEOPOT. DEPENDING UPON Tv=tpnv() C**** phnhat() = (full level) GEOPOT. DEPENDING UPON T"hat" C**** VERTICAL COORDINATE (HYBRID/SIGMA) FUNCTIONS AND CONSTANTS if(hybrid)then fhyb=1.0 do 2076 k=1,nl do 2076 mg=1,ln2 fff1(mg,k)=muf(mg,k)/prf(mg,k)*sig(k) fff2(mg,k)=muf(mg,k)/prf(mg,k)*bnf(k) 2076 fff3(mg,k)=prf(mg,k)/muf(mg,k) else fhyb=0.0 do 2078 k=1,nl do 2078 mg=1,ln2 fff1(mg,k)=1.0 fff2(mg,k)=1.0 2078 fff3(mg,k)=sig(k) endif do 2080 k=1,nl do 2080 mg=1,ln2 txn1(mg,k)=fchen*fff1(mg,k)*tdn(mg,k) & +(1.0-fchen)*(fff1(mg,k)*tpn(mg,k)-tmean(k)) txn2(mg,k)=fff1(mg,k)*(tpnv(mg,k)-tpn(mg,k)) 2080 txn3(mg,k)=fff1(mg,k)*tpnv(mg,k) c**** USE MATRIX TO EXPRESS GEOPOT.HEIGHT AS TEMPERATURE do 207 mg=1,ln2 207 phnh(mg,1)=0.0 do 122 k=1,nl do 208 mg=1,ln2 phnv(mg,k) = 0.0 phn(mg,k)=phist(mg) phnh(mg,k+1)=0.0 208 phnhat(mg,k)=0.0 do 122 j=1,nl do 209 mg=1,ln2 phnv(mg,k)=phnv(mg,k)+am(k,j)*txn2(mg,j) phn(mg,k)=phn(mg,k)+am(k,j)*txn1(mg,j) phnh(mg,k+1)=phnh(mg,k+1)+amh(k,j)*txn3(mg,j) 209 phnhat(mg,k)=phnhat(mg,k)+am(k,j)*ten(mg,j) 122 continue if(lstat)then c Full level GPH computed for Stats purposes only do mg=1,ln2 phisf(mg)=max(0.0,grav*z4(mg)) enddo do k=1,nl do mg=1,ln2 phnf(mg,k)=phisf(mg) enddo do j=1,nl do mg=1,ln2 phnf(mg,k)=phnf(mg,k)+am(k,j)*txn3(mg,j) enddo enddo enddo endif C**** COMPUTE NON-LINEAR FLUX TERMS AT EACH LEVEL C** SET UP VERTICAL FLUX TERMS do 214 mg=1,ln2 temu(mg,1)=-uh(mg,2)/dsk(1) temv(mg,1)=-vh(mg,2)/dsk(1) temt(mg,1)=-th(mg,2)/dsk(1) temu(mg,nl)=uh(mg,nl)/dsk(nl) temv(mg,nl)=vh(mg,nl)/dsk(nl) 214 temt(mg,nl)=th(mg,nl)/dsk(nl) do 128 k=2,nlm do 128 mg=1,ln2 temu(mg,k)=(uh(mg,k)-uh(mg,k+1))/dsk(k) temv(mg,k)=(vh(mg,k)-vh(mg,k+1))/dsk(k) 128 temt(mg,k)=(th(mg,k)-th(mg,k+1))/dsk(k) C* C** ADD ON VARIOUS COMPONENTS OF NONLINEAR TERMS do mg=1,lon ns=1 tomgns(mg )=tomg*coa(lg)*(3.0-2.0*ns) ns=2 tomgns(mg+lon)=tomg*coa(lg)*(3.0-2.0*ns) enddo do 1302 k=1,nl do 1302 mg=1,ln2 fon(mg,k)=un(mg,k)*tdn(mg,k) gon(mg,k)=vn(mg,k)*tdn(mg,k) upnk=upn(mg,k) vpnk=vpn(mg,k) vonk=von(mg,k)+dbdnf(k)*(upnk*cpn(mg)-vpnk*pln(mg)) & *csqlg+muf(mg,k)*tomgns(mg) etem=(upnk*upnk+vpnk*vpnk)*0.5*csqlg-emean(k) C**** USE VIRTUAL TEMPERATURE AND CHANGE FOR "CHENG" temz1=fchen*(fff2(mg,k)*(tpnv(mg,k)+(tdn(mg,k)-tpn(mg,k))) & -t00ec) & +(1.0-fchen)*(fff2(mg,k)*tpnv(mg,k)-tmean(k)) tem=( rdry*temz1-dbdnf(k)*(phn(mg,k)+phnv(mg,k)))/eradsq & -etem*dbdnf(k) aon(mg,k)=dvn(mg,k)*vpn(mg,k)-son(mg,k) & +vonk*upnk+temv(mg,k)+tem*cpn(mg) bon(mg,k)=ron(mg,k)-dvn(mg,k)*upn(mg,k) & +vonk*vpnk-temu(mg,k)-tem*pln(mg) C**** ADD THE phnv IN HERE eon(mg,k)=(etem+(phnv(mg,k)+phist(mg))/eradsq)*muf(mg,k) & +fhyb*(muf(mg,k)*(phn(mg,k)-phist(mg))-phnhat(mg,k))/eradsq delt2=(upnk*pln(mg)+vpnk*cpn(mg))*csqlg C DVGG(MG,LG,INS)=DVGG(MG,LG,INS) C * +DSKD(K)*(DVN(MG,K)-DELT2)/PN(MG) C**** dpsdt FOR jmcgslt ONLY dpsdtk(mg,k)=delt2-dvhn(mg,1) c wsx=delt2-(dvfn(mg,k)/sig(k)) !omega/sigma wsx=fff2(mg,k)*delt2-(dvfn(mg,k)/fff3(mg,k)) !mu*omega/p c omega(mg,k,lg)=bnf(k)*delt2 - dvfn(mg,k) !mb/s c write(27,'(2f13.5)')float(k),omega C**** CHANGE FOR "CHENG" hon(mg,k)=cappa*(fff2(mg,k)*tpnv(mg,k)*delt2+dvfn(mg,k)*( & (fchen*t00ec+(1.0-fchen)*tmean(k))/sig(k) & +(phnh(mg,k)-phnh(mg,k+1))/ & (rdry*dsk(k))) +fchen*wsx*(tdn(mg,k)-tpn(mg,k)+t0ec)) & -temt(mg,k) 1302 continue do 1303 k=1,nl do 1303 mg=1,lon sdot(mg,k,lgn)=0.5*(wbn(mg,k)+wbn(mg,k+1))/muf(mg,k) sdot(mg,k,lg )=0.5*(wbn(mg+lon,k)+wbn(mg+lon,k+1))/muf(mg+lon,k) dpsdt(mg,lgn,k)=dpsdtk(mg,k) dpsdt(mg,lg ,k)=dpsdtk(mg+lon,k) c sdot2(mg,k,lgn)=sdot(mg,k,lgn)+omegac(mg,k,lg)/muf(mg,k) !Approx c sdot2(mg,k,lg)=sdot(mg,k,lg)+omegac(mg+lon,k,lg)/muf(mg+lon,k) 1303 continue C**** GATHER STATISTICS EVERY 1/4 DAY ONLY zonpsl(1)=0.0 zonpsl(2)=0.0 if(.not.lstat)return flong=1.0/lon c=grav/6.5e-3 conr=c/rdry C* CORRECT PST TO MEAN SEA LEVEL PSL if(hybrid)then do 2192 mg=1,ln2 pxlev=pn(mg)-90.0 k=1 999 k=k+1 if(pxlev.lt.prf(mg,k))go to 999 plk(mg,2)=prf(mg,k) plk(mg,1)=prf(mg,k-1) tlk(mg,2)=tpn(mg,k) 2192 tlk(mg,1)=tpn(mg,k-1) do 2194 mg=1,ln2 pxlev=pn(mg)-90.0 txlev=(tlk(mg,2)*alog(pxlev/plk(mg,1)) + & tlk(mg,1)*alog(plk(mg,2)/pxlev))/alog(plk(mg,2)/plk(mg,1)) con=(pxlev/pn(mg))**(rdry/c) /c psl(mg)=pn(mg)*(1.+con*grav*z4(mg)/txlev) **conr 2194 dvgg(mg,lg)=psl(mg) else con=sig(2)**(rdry/c) /c do 219 mg=1,ln2 psl(mg)=pn(mg)*(1.+con*grav*z4(mg)/tpn(mg,2)) **conr 219 dvgg(mg,lg)=psl(mg) endif C**** C**** COLLECT DYNAMICAL STATS FOR DAILY PRINT OUT. C*** COLLECT MERIDIONAL TRANSPORT STATS (NORTHWARD FLUXES). C ZONALLY AVERAGED NOTHWARD FLUXES COMPUTED AT EACH LEVEL C EG : ZONAL AVG (V.T) AND ZONAL AVG (V). ZONAL AVG (T) C THIS WILL YIELD MEAN + TRANSIENT FLUXES. C DONE FOR V.T, V.U, V.Q C ALSO COLLECTING ZONAL AVG (P*.V) FOR STREAM FN. do 4444 ns=1,2 ma=(ns-1)*lon do 631 k=1,nl uz=0.0 vz=0.0 tz=0.0 pvtbz=0.0 pvubz=0.0 psvz=0.0 do 221 mg=1,lon uz=uz+upn(mg+ma,k) vz=vz+vpn(mg+ma,k) tz=tz+tpn(mg+ma,k) pvtbz=pvtbz+vpn(mg+ma,k)*tpn(mg+ma,k) pvubz=pvubz+vpn(mg+ma,k)*upn(mg+ma,k) 221 psvz=psvz+vn(mg+ma,k) uz=uz*siai*flong vz=vz*siai*flong tz=tz*flong u(lg,ns,k)=u(lg,ns,k)+uz v(lg,ns,k)=v(lg,ns,k)+vz t(lg,ns,k)=t(lg,ns,k)+tz pvbtb(lg,ns,k)=pvbtb(lg,ns,k)+vz*tz pvbub(lg,ns,k)=pvbub(lg,ns,k)+vz*uz pvtb(lg,ns,k)=pvtb(lg,ns,k)+pvtbz*siai*flong pvub(lg,ns,k)=pvub(lg,ns,k)+pvubz*siai*siai*flong psv(lg,ns,k)=psv(lg,ns,k)+psvz*siai*flong if(k.gt.nl)go to 631 qz=0.0 pvqbz=0.0 do 223 mg=1,lon qz=qz+rpn(mg+ma,k) 223 pvqbz=pvqbz+vpn(mg+ma,k)*rpn(mg+ma,k) pvbqb(lg,ns,k)=pvbqb(lg,ns,k)+vz*(qz*flong) pvqb(lg,ns,k)=pvqb(lg,ns,k)+pvqbz*siai*flong 631 continue do 632 k=1,nlm wzk=0.0 do 224 mg=1,lon 224 wzk=wzk+wbn(mg+ma,k+1) wzk=wzk*flong C ZONAL VERT VELY P*.SIGMA(DOT) pssd(lg,ns,k)=pssd(lg,ns,k)+wzk 632 ww(lg,ns,k)=ww(lg,ns,k)+wzk psz=0.0 do 225 mg=1,lon 225 psz=psz+psl(mg+ma) zonpsl(ns)=psz ps(lg,ns)=ps(lg,ns)+psz*flong 4444 continue return end