!
! $Id: pdyn0.F90,v 1.1 2005/10/18 20:10:16 mbhat Exp $
!
!---(pdiag.f)-------generic CTM shell from UCIrvine (p-code 1.0, 1/31/95)
!---subroutines: DYN0, GRDSET, AIRSET, PFILTR
SUBROUTINE DYN0(LOOP) 1,1
!---------QDYN0 sets up the advective fields
!---CURRENTLY THIS IS CALLED FROM MAIN WITH
!--- LOOP=0 = SETUP, OR LOOP=1 = CONTINUE
!---PLAN TO ELIMINATE LOOP=0 CALL!
!---THIS SET OF ARRAYS & COMMONS IS A TEMPORARY FIX TO MAKE THE
!--- NEW VERSION OF -DYN0- BASED ON ECMWF GRID, COMPATIBLE WITH
!--- THE OLD SET OF COMMONS, ETC TO BE USED IN FIRST LLNL VERSION (QZ=2.0)
INTEGER LOOP
LOGICAL LSP,LNP,LEW
REAL*8 SUMA,SUMP,SUMQ,SUMU,SUMV,SUMW
REAL*8 AIRD,AIRNEW,AIRX, XYZA,XYZB,XYA,XYB
REAL*8 SUMXYZ,SUMAD0,SUMAW0,AIRWET,AIRH2O,AIRDRY,AIRQKG
REAL*8 SUMAQ(36,24)
REAL*8 PERR(36,24),MERR(36,24),AX(37,24),BX(36,25),DTWIND
!---
COMMON /CCGRID/ &
& AIRD(36,24,21),AIRNEW(36,24,21),AIRX(36,24,21), &
& XYZA(36,24,21),XYZB(36,24,21),XYA(36,24),XYB(36,24), &
& AIRDRY,AIRQKG, &
& YGRD(24),YDGRD(24),YEDG(25),YDEDG(25),XGRD(36), &
& XDGRD(36),XEDG(37),XDEDG(37),DLONG,AREAXY(36,24), &
& DISTX(25),DISTY(24), ETAA(22),ETAB(22), &
& UUU(37,24,21),VVV(36,25,21),PS(36,24), &
& AIRQG(36,24,21),AIRQ(36,24,21),PCTM(36,24), IMEPZ(24)
!---
!----------------------------------------
# include "uci.h"
!xx$$
save
!xx$$
!-----------------------------------------------------------------------
!----------------- UNITS OF AIR MASS AND TRACER = (kg) ----------------
!----Air mass (kg) is given by area (m^2) * pressure thickness (Pa) / g0
!---- AREAXY(I,J) = area of (I,J) (m^2)
!---- PS(I,J) = surf pressure (Pa) as reported by GCM
!---- AIRQ(I,J,K) = specific humidity of grid box (kg H2O / kg wet air) per GCM
!---- AIRQKG = local (I,J) kg of H2O at each level based on GCM P's & q's
!---- = PS(I,J)) * AIRQ(I,J,K)
!---- AIRD(I,J,K) = dry-air mass (kg) in each box as calculated in CTM
!---- at the beginning of each time step, updated at end of DYN0.
!---- PCTM(I,J) = inferred wet-air (total) surf press (Pa) calc. in CTM
!---- (using SUMAQ & AIRD-X-NEW)
!---- AIRNEW(I,J,K) = new dry-air mass in each CTM box after horizontal
!---- divergence (ALFA+BETA) over time step DTWIND (sec)
!---- AIRX(I,J,K) = expected dry-air mass in each CTM box after calculating
!---- the vertical divergence (GAMA) (also used for GCM dry mass)
!---- = XYZA(I,J,K) + XYZB(I,J,K)*PCTM(I,J) - AIRQKG
!----
!----Assume that we have "wet-air" mass fluxes across each boundary
!---- U(I,J,K) ==> (I,J,K) ==> U(I+1,J,K) (kg/s)
!---- V(I,J,K) ==> (I,J,K) ==> V(I,J+1,K) (kg/s)
!--- DTWIND = time step (sec) that applies to the averaged wind fields
!--- i.e., the time between successive pressures.
!----Convert to "dry-air" mass flux in/out of box using average AIRQ at boundary
!---- ALFA(I,J,K) ==> (I,J,K) ==> ALFA(I+1,J,K) (kg/s)
!---- BETA(I,J,K) ==> (I,J,K) ==> BETA(I,J+1,K) (kg/s)
!----
!----Calculate convergence in each layer of dry air, compare with expected
!---- dry air mass (AIRX) and then calculate vertical dry-mass fluxes
!---- GAMA(I,J,K) ==> (I,J,K) ==> GAMA(I,J,K+1) (kg/s)
!----
!----Horizontal pressure filter adjusts U & V to reduce error in (PCTM - PS)
!---- U + pressure filter ==> U#, V + filter ==> V# (temporary)
!---- The pressure filter does nearest neighbor flux (adjusting ALFA/BETA)
!----
!----Time step for each of the operator-splitting calculations in U/V/W
!---- is determined by the number of U-V-W calls per time step (NDYN)
!---- and DTALFA/DTBETA/DTGAMA must be set for DYN2U/V/W (in main driver)
!-----------------------------------------------------------------------
!----
!----Note that K->K+1 is downward (increasing pressure) and that boundaries:
!---- GAMA(I,J,1) = GAMA(I,J,KM+1) = 0 no flux across upper/lower boundaries
!---- BETA(I,1,K) = BETA(I,JM+1,K) = 0 no flux at S & N poles
!---- ALFA(1,J,K) = ALFA(IM+1,J,K) is NOT ZERO, but cyclic
!----Dimensions for ALFA, BETA, GAMA are extended by +1 beyond grid to allow
!---- simple formulation of fluxes in/out of final grid box.
!
!-----GCM input U,V,AIRQG,PSG is ALWAYS of GLOBAL dimensions (IMG x JMG x KMG)
!-----Indices of ALFA, BETA, GAMA, AIRQ & PS are always LOCAL (IM x JM x KM)
!-----Indices of tracer (STT), and diagnostics are local (w.r.t. WINDOW)
!-----WINDOW calculations are defined by an offset and size
!----- I0 .ge.0 and IM+I0 .le. IMG
!----- J0 .ge.0 and JM+J0 .le. JMG
!----- K0 .ge.0 and KM+K0 .le. KMG
!-----The WINDOW calculation must allow for a boundary layer of grid boxes
!----- IG(abs. coords) = IW(in window) + I0
!----- JG(abs. coords) = JW(in window) + J0
!----- KG(abs. coords) = KW(in window) + K0
!----- vertical window (NEW) allows for an upper boundary with flow across
!----- it and specified mixing ratio b.c.'s at KG = K0
!-------------------------------------------------------------------------
!
KWACC = KWACC+1
!---------NREAD = no. hours between pressure fields (U & V are averaging time)
DTWIND = (3600*NREAD)
G0 = 9.81D0
IMG = IMX
JMG = JMX
KM = LM
!---------First need to map the full-grid values (__G) onto local(WINDOW) arrays
!---------NB need to have IMG=IM, JMG=JM, KMG=KM for full global calc.
DO J=1,JM
DO I=1,IM
II0 = MOD(I+I0-1, IMG) + 1
!---------original for q-code: PS(I,J) = PSG(II0,J+J0)--------------
PS(I,J) = P(II0,J+J0)
DO K=1,KM
AIRQ(I,J,K) = AIRQG(II0,J+J0,K)
ENDDO
ENDDO
ENDDO
!---------TEMP FIX OF WINDS: NO WINDOWS ALLOWED
G100 = 100.0/G0
DO L=1,LM
DO J=1,JMG
DO I=1,IMG
UUU(I,J,L) = U(I,J,L)*G100
ENDDO
UUU(IMG+1,J,L) = UUU(1,J,L)
ENDDO
ENDDO
!---------TEMP FIX OF WINDS: NO WINDOWS ALLOWED
DO L=1,LM
DO I=1,IMG
DO J=2,JMG
VVV(I,J,L) = V(I,J,L)*G100
ENDDO
VVV(I,1,L) = 0.D0
VVV(I,JMG+1,L) = 0.D0
ENDDO
ENDDO
!---------PROBLEMS WITH POLAR NOISE IN GISS 21-LAYER: temporary
!---V-8.2----scale V-field at polar box, proportional to area J=1/(1+2)
!---------divide polar convergence over J=1&2 and J=JMX-1&JMX (assume symmetry)
DO L=1,LM
DO I=1,IMG
VVV(I,2,L) = VVV(I,3,L)*AREAXY(I,1)/(AREAXY(I,1)+AREAXY(I,2))
VVV(I,JMG,L) = VVV(I,JMG-1,L)* &
& AREAXY(I,JMG)/(AREAXY(I,JMG)+AREAXY(I,JMG-1))
ENDDO
ENDDO
!--------------------------END OF TEMP FIX------------------------------
DO K=1,KM
DO J=1,JM
DO I=1,IM+1
II0 = MOD(I+I0-1, IMG) + 1
ALFA(I,J,K) = UUU(II0,J+J0,K)
ENDDO
ENDDO
DO I=1,IM
II0 = MOD(I+I0-1, IMG) + 1
DO J=1,JM+1
BETA(I,J,K) = VVV(II0,J+J0,K)
ENDDO
ENDDO
ENDDO
!---------DRY-AIR: re-define ALFA & BETA to be dry-air flux (kg/s),
!---------cannot go beyond boundaries of WINDOW box:
DO K=1,KM
DO J=1,JM
IF(J+J0.EQ.1 .OR. J+J0.EQ.JMG) THEN
!---------safety zero of ALFAs in polar boxes
DO I=1,IM+1
ALFA(I,J,K) = 0.0
ENDDO
ELSE
DO I=2,IM
AIRQAV = 0.5*(AIRQ(I-1,J,K)+AIRQ(I,J,K))
ALFA(I,J,K) = ALFA(I,J,K)*(1.-AIRQAV)
ENDDO
IF(IM.LT.IMG) THEN ! WINDOW
ALFA(1,J,K) = ALFA(1,J,K)*(1.-AIRQ(1,J,K))
ALFA(IM+1,J,K) = ALFA(IM+1,J,K)*(1.-AIRQ(IM,J,K))
ELSE ! full cyclic grid
AIRQAV = 0.5*(AIRQ(1,J,K)+AIRQ(IM,J,K))
ALFA(1,J,K) = ALFA(1,J,K)*(1.-AIRQAV)
ALFA(IM+1,J,K) = ALFA(1,J,K)
ENDIF
ENDIF
ENDDO ! end of J loop
DO J=2,JM ! interior J points
DO I=1,IM
AIRQAV = 0.5*(AIRQ(I,J-1,K)+AIRQ(I,J,K))
BETA(I,J,K) = BETA(I,J,K)*(1.-AIRQAV)
ENDDO
ENDDO ! end of J loop
IF(J0.EQ.0) THEN !flux into S.Pole = 0
DO I=1,IM
BETA(I,1,K) = 0.0
ENDDO
ELSE
DO I=1,IM !flux across southernmost WINDOW boundary
BETA(I,1,K) = BETA(I,1,K)*(1.-AIRQ(I,1,K))
ENDDO
ENDIF
IF(JM+J0.EQ.JMG) THEN !flux out of N.Pole = 0
DO I=1,IM
BETA(I,JM+1,K) = 0.0
ENDDO
ELSE
DO I=1,IM
BETA(I,JM+1,K) = BETA(I,JM+1,K)*(1.-AIRQ(I,JM,K))
ENDDO
ENDIF
ENDDO ! end of K loop
!----EPZ'S:
!------Average GCM quantities: PS, AIRQ, BETA
DO J=1,JM
IMZ = IMEPZ(J+J0)
IF(IMZ.GT.1) THEN
!----average PS over EPz's
DO IEPZ=1,IM,IMZ ! # EPZ's = IM/IMZ
SUMP = 0.D0
DO I=IEPZ,IEPZ+IMZ-1 ! loop over I's in each EPZ
SUMP = SUMP + PS(I,J)
ENDDO
SUMP = SUMP/REAL(IMZ) !SUMP=pressure (w-air) in box(I,J)
DO I=IEPZ,IEPZ+IMZ-1
PS(I,J) = SUMP
ENDDO
ENDDO
!----average AIRQ over EPZ's
DO K=1,KM
DO IEPZ=1,IM,IMZ ! # EPZ's = IM/IMZ
SUMA = 0.D0
SUMQ = 0.D0
DO I=IEPZ,IEPZ+IMZ-1
!---- have not bothered to offset I since XYZA/B independent of I
AWET = XYZA(I,J,K) + PS(I,J)*XYZB(I,J,K)
SUMA = SUMA + AWET
SUMQ = SUMQ + AIRQ(I,J,K)*AWET
ENDDO
SUMQ = SUMQ/SUMA
DO I=IEPZ,IEPZ+IMZ-1
AIRQ(I,J,K) = SUMQ
ENDDO
ENDDO
ENDDO ! end of K loop
ENDIF ! check on IMZ
ENDDO ! end of J loop
!----Diagnostics and use of PW() for z-code overlap
DO J=1,JM
DO I=1,IM
PW(I,J) = PS(I,J) - PTOP
AIJ(I,J,1) = AIJ(I,J,1) + PS(I,J)
ENDDO
ENDDO
!---------SUMAQ(I,J): column integral of water (kg)---------------------
!---------check on air mass---------------------------------------------
DO J=1,JM
DO I=1,IM
SUMAQ(I,J) = 0.D0
ENDDO
ENDDO
SUMAD0 = 0.D0
SUMAW0 = 0.D0
DO K=1,KM
DO J=1,JM
DO I=1,IM
AIRWET = XYZA(I,J,K) + PS(I,J)*XYZB(I,J,K)
AIRH2O = AIRQ(I,J,K)*AIRWET
SUMAQ(I,J) = SUMAQ(I,J) + AIRH2O
SUMAD0 = SUMAD0 + AIRWET
SUMAW0 = SUMAW0 + AIRH2O
ENDDO
ENDDO
ENDDO
SUMAD0 = SUMAD0 - SUMAW0
IF(LPRT) WRITE(lun13,'(A,1P,2E12.5,A,E12.5)') &
& ' ---dry air: init/now',AIRDRY,SUMAD0,' ---H2O: now',SUMAW0
!---------smooth BETA's over EPZ's: (use closest to equator value for EPZ)
!---------i.e., BETA(3) is between boxes (2) and (3), average over IMEPZ(3)
DO J=1,JM+1
IF(2*(J+J0).LE.JMG) THEN
IMZ = IMEPZ(J+J0) ! S.Hem.
ELSE
IMZ = IMEPZ(J+J0-1) ! N.Hem.
ENDIF
IF(IMZ.GT.1 .AND. .NOT.(J+J0.EQ.1 .OR. J+J0.EQ.JMG+1)) THEN !skip poles
DO K=1,KM
DO IEPZ=1,IM,IMZ ! # EPZ's = IM/IMZ
SUMV = 0.D0
DO I=IEPZ,IEPZ+IMZ-1
SUMV = SUMV + BETA(I,J,K)
ENDDO
SUMV = SUMV/REAL(IMZ)
DO I=IEPZ,IEPZ+IMZ-1
BETA(I,J,K) = SUMV
ENDDO
ENDDO
ENDDO
ENDIF ! check on IMZ
ENDDO ! end of J loop
!---CTM CONVERGENCE: DRY-AIR mass convergence into each CTM grid box:
DO K=1,KM
DO J=1,JM
DO I=1,IM
!-------*****AD() IS THE ONLY AIR MASS CALC IN Z-CODE
AIRD(I,J,K) = AD(I,J,K)
AIRNEW(I,J,K) = AIRD(I,J,K) + DTWIND* &
& (ALFA(I,J,K)-ALFA(I+1,J,K)+BETA(I,J,K)-BETA(I,J+1,K))
ENDDO
ENDDO
ENDDO
!---reset ALFAs at interior pts of EPZ's to guarantee equal masses:
DO K=1,KM
DO J=1,JM
IF(J+J0.NE.1 .AND. J+J0.NE.JMG) THEN !not for poles
IMZ = IMEPZ(J+J0)
IF(IMZ.GT.1) THEN
DO IEPZ=1,IM,IMZ !IEPZ strides thru EPZs: 1, 1+IMZ, 1+2*IMZ,..
SUMA = 0.D0
DO I=IEPZ,IEPZ+IMZ-1 !I all boxes in a single EPZ
SUMA = SUMA + AIRNEW(I,J,K)
ENDDO
SUMA = SUMA/REAL(IMZ)
ALFAX = (AIRNEW(IEPZ,J,K)-SUMA)/DTWIND
DO I=IEPZ+1,IEPZ+IMZ-1
ALFA(I,J,K) = ALFA(I,J,K) + ALFAX
ALFAX = ALFAX + (AIRNEW(I,J,K) - SUMA)/DTWIND
ENDDO
DO I=IEPZ,IEPZ+IMZ-1
AIRNEW(I,J,K) = SUMA
ENDDO
ENDDO
ENDIF
ELSE ! poles, just average AIRNEW
SUMA = 0.D0
DO I=1,IM
SUMA = SUMA + AIRNEW(I,J,K)
ENDDO
SUMA = SUMA/REAL(IM)
DO I=1,IM
AIRNEW(I,J,K) = SUMA
ENDDO
ENDIF
ENDDO !J
ENDDO !K
!-----------------------------------------------------------------------
!--------------------------BEGIN FILTER of PRESSURE ERRORS--------------
!-----------------------------------------------------------------------
!---Define the error in surface pressure PERR expected at end of time step
!------filter by error in adjacent boxes, weight by areas, adjust ALFA & BETA
!
IF(.NOT.LPFILT) GOTO 101
!
!---PCTM(I,J)= new CTM wet-air column based on dry-air convergence (Pascals)
!---PERR(I,J)= pressure-error between CTM-GCM at new time (before filter)
DO J=1,JM
DO I=1,IM
SUMA = 0.D0
DO K=1,KM
SUMA = SUMA + AIRNEW(I,J,K)
ENDDO
PCTM(I,J) = (SUMA + SUMAQ(I,J) - XYA(I,J))/XYB(I,J)
PERR(I,J) = PCTM(I,J) - PS(I,J)
MERR(I,J) = PERR(I,J)*AREAXY(I,J+J0)*100.E0/G0
ENDDO
ENDDO
!-----------------------------------------------------------------------
LSP = J0.EQ.0
LNP = JM+J0.EQ.JMG
LEW = IM.EQ.IMG
! CALL PFILTR(MERR,AX,BX,IM,JM,IMG,JMG,IMEPZ(J0+1),LSP,LNP,LEW)
CALL PFILTR
(MERR,AX,BX,IM,JM,IMG,JMG,IMEPZ,LSP,LNP,LEW)
!---------Calculate corrections to ALFA & BETA from AX and BX:----------
DO J=1,JM
IF(J+J0.GT.1 .AND. J+J0.LT.JMG) THEN ! NOT Poles
DO I=1,IM+1
IIX = MIN(I,IM)
UFILT = AX(I,J)/(XYB(IIX,J)*DTWIND)
DO K=1,KM
ALFA(I,J,K) = ALFA(I,J,K) + UFILT*XYZB(IIX,J,K)
ENDDO
ENDDO
ENDIF
ENDDO
DO J=1,JM+1
JJX = J
IF(J+J.GT.JM) JJX=J-1
DO I=1,IM
VFILT = BX(I,J)/(XYB(I,JJX)*DTWIND)
DO K=1,KM
BETA(I,J,K) = BETA(I,J,K) + VFILT*XYZB(I,JJX,K)
ENDDO
ENDDO
ENDDO
!----Calculate the corrected AIRNEW's & PCTM after P-filter:
!---- has changed ALFA+BETAs and ctm surface pressure (PCTM)
DO J=1,JM
DO I=1,IM
DO K=1,KM
AIRNEW(I,J,K) = AIRD(I,J,K) + DTWIND* &
& (ALFA(I,J,K)-ALFA(I+1,J,K)+BETA(I,J,K)-BETA(I,J+1,K))
ENDDO
ENDDO
ENDDO
!---reset ALFAs at interior pts of EPZ's to guarantee equal masses:
DO K=1,KM
DO J=1,JM
IF(J+J0.NE.1 .AND. J+J0.NE.JMG) THEN !not for poles
IMZ = IMEPZ(J+J0)
IF(IMZ.GT.1) THEN
DO IEPZ=1,IM,IMZ !IEPZ strides thru EPZs: 1, 1+IMZ, 1+2*IMZ,..
SUMA = 0.D0
DO I=IEPZ,IEPZ+IMZ-1 !I all boxes in a single EPZ
SUMA = SUMA + AIRNEW(I,J,K)
ENDDO
SUMA = SUMA/REAL(IMZ)
ALFAX = (AIRNEW(IEPZ,J,K)-SUMA)/DTWIND
DO I=IEPZ+1,IEPZ+IMZ-1
ALFA(I,J,K) = ALFA(I,J,K) + ALFAX
ALFAX = ALFAX + (AIRNEW(I,J,K) - SUMA)/DTWIND
ENDDO
DO I=IEPZ,IEPZ+IMZ-1
AIRNEW(I,J,K) = SUMA
ENDDO
ENDDO
ENDIF
ELSE ! poles, just average AIRNEW
SUMA = 0.D0
DO I=1,IM
SUMA = SUMA + AIRNEW(I,J,K)
ENDDO
SUMA = SUMA/REAL(IM)
DO I=1,IM
AIRNEW(I,J,K) = SUMA
ENDDO
ENDIF
ENDDO !J
ENDDO !K
!--------------------------end of pressure filter-----------------------
101 CONTINUE
!---GAMA: redistribute the new dry-air mass consistent with the
!--- new CTM surface pressure, rigid upper b.c., no change in PCTM
!---AIRX(I,J,K) = dry-air mass expected, based on PCTM
!---PCTM(I,J) & PERR(I,J)
DO J=1,JM
DO I=1,IM
SUMA = 0.D0
DO K=1,KM
SUMA = SUMA + AIRNEW(I,J,K)
ENDDO
PCTM8 = (SUMA + SUMAQ(I,J) - XYA(I,J))/XYB(I,J)
PCTM(I,J) = PCTM8
PERR(I,J) = PCTM(I,J) - PS(I,J)
DO K=1,KM
AIRQKG = AIRQ(I,J,K)*(XYZA(I,J,K)+XYZB(I,J,K)*PS(I,J))
AIRX(I,J,K) = XYZA(I,J,K) + PCTM8*XYZB(I,J,K) - AIRQKG
ENDDO
ENDDO
ENDDO
!---GAMA from top down to be consistent with AIRX, AIRNEW not reset!
DO J=1,JM
DO I=1,IM
GAMA(I,J,KM+1) = 0.D0
DO K=KM,2,-1
GAMA(I,J,K) = GAMA(I,J,K+1) - (AIRNEW(I,J,K) - AIRX(I,J,K))
ENDDO
!---N.B. GAMA(I,J,1) will not be exactly ZERO, but it must be set so!
GAMA(I,J,1) = 0.D0
!---REPLACE SECTION for K=TOP:BOT for EC-model
! GAMA(I,J,1) = 0.E0
! DO K=1,KM-1
! GAMA(I,J,K+1) = GAMA(I,J,K) + (AIRNEW(I,J,K) - AIRX(I,J,K))
! ENDDO
! GAMA(I,J,KM+1) = 0.E0
!---scale GAMA from kg over the wind-field time interval to kg/sec
DO K=2,KM
GAMA(I,J,K) = GAMA(I,J,K)/DTWIND
ENDDO
ENDDO
ENDDO
! WRITE OUT ENTIRE ARRAY
IF(.FALSE.) THEN
write(6,*) 'writing out pctm array'
write(16,'(e16.6)') ((pctm(i,j),i=1,im),j=1,jm)
write(6,*) 'writing out alfa array'
write(16,'(e16.6)') (((alfa(i,j,k),i=1,im),j=1,jm ),k=1,km )
write(6,*) 'writing out beta array'
write(16,'(e16.6)') (((beta(i,j,k),i=1,im),j=1,jm ),k=1,km )
write(6,*) 'writing out gama array'
write(16,'(e16.6)') (((gama(i,j,k),i=1,im),j=1,jm ),k=1,km )
call Esm_Flush(16)
ENDIF
!---DIAGNOSTICS:
!--- mean winds N.B. diagnostics are DRY-air winds ONLY!
!
!----accumulate zonal average U (m/s)
DO K=1,KM
DO J=1,JM
SUMU = 0.D0
DO I=1,IM
SUMU = SUMU + ALFA(I,J,K)/AIRD(I,J,K)
ENDDO
AJL(J,K,1) = AJL(J,K,1) + SUMU*DISTX(J+J0)/REAL(IM)
ENDDO
ENDDO
!
!----accumulate zonal average V (m/s)
DO K=1,KM
DO J=2,JM
SUMV = 0.D0
DO I=1,IM
SUMV = SUMV + BETA(I,J,K)/(AIRD(I,J-1,K)+AIRD(I,J,K))
ENDDO
AJL(J,K,2) = AJL(J,K,2) + SUMV*(DISTY(J+J0-1) + &
& DISTY(J+J0))/REAL(IM)
ENDDO ! end of loop J=2,JM, now do J=1 & J=JM+1
SUMV = 0.0
DO I=1,IM
SUMV = SUMV + BETA(I,1,K)/AIRD(I,1,K)
ENDDO
AJL(1,K,2)=AJL(1,K,2)+SUMV*DISTY(1+J0)/REAL(IM)
IF(JM+1.LT.JMG) THEN
SUMV = 0.D0
DO I=1,IM
SUMV = SUMV + BETA(I,JM+1,K)/AIRD(I,JM,K)
ENDDO
AJL(JM+1,K,2)=AJL(JM+1,K,2)+SUMV*DISTY(JM+J0)/REAL(IM)
ENDIF
ENDDO
!
!----accumulate zonal average Omega (Pascal/sec)
DO K=1,KM
! WRITE(lun13,'(12HNTAU/K/PLEV=,I5,I2,1PE12.5)') NTAU,K,PLEV(K,1)
DO J=1,JM
SUMA = G0 /(100. *PLEV(K,1) *AREAXY(1,J+J0) *REAL(IM))
SUMW = 0.D0
DO I=1,IM
SUMW = SUMW + GAMA(I,J,K)
ENDDO
AJL(J,K,3) = AJL(J,K,3) + SUMW *SUMA
ENDDO
ENDDO
!
RETURN
END
!
!---------subroutine GRDSET--------------------------------------------
!
SUBROUTINE GRDSET 1
!
!---------GRDSET sets up the grids
!---------THIS SET OF ARRAYS & COMMONS IS A TEMPORARY FIX TO MAKE THE
!---------NEW VERSION OF -DYN0- BASED ON ECMWF GRID, COMPATIBLE WITH
!---------THE OLD SET OF COMMONS, ETC TO BE USED IN FIRST LLNL VERSION (QZ=2.0)
REAL*8 AIRD,AIRNEW,AIRX, XYZA,XYZB,XYA,XYB
REAL*8 SUMXYZ,SUMAD0,SUMAW0,AIRWET,AIRH2O,AIRDRY,AIRQKG
!
COMMON /CCGRID/ &
& AIRD(36,24,21),AIRNEW(36,24,21),AIRX(36,24,21), &
& XYZA(36,24,21),XYZB(36,24,21),XYA(36,24),XYB(36,24), &
& AIRDRY,AIRQKG, &
& YGRD(24),YDGRD(24),YEDG(25),YDEDG(25),XGRD(36), &
& XDGRD(36),XEDG(37),XDEDG(37),DLONG,AREAXY(36,24), &
& DISTX(25),DISTY(24), ETAA(22),ETAB(22), &
& UUU(37,24,21),VVV(36,25,21),PS(36,24), &
& AIRQG(36,24,21),AIRQ(36,24,21),PCTM(36,24), IMEPZ(24)
!-----------------------------------------------------------------------
# include "uci.h"
!xx$$
save
!xx$$
!-----------------------------------------------------------------------
!---------AREAXY(I,J) = area of (I,J) (m^2)
!---------pressure at I,J,K is ETTA(K) + ETTB(K)*PS(I,J)---------------
!-------------------------------------------------------------------------
!
!---------this initialization is temporary, transition to new notation
IMG = IMX
JMG = JMX
KM = LM
DO J=1,JMG
IMEPZ(J) = IMRSLV(J)
ENDDO
!----------------------GRID---------------------------------------------
!---------fix for GISS 21 layer
!---------ECMWF: A0 = 6371000.d0
!---------ECMWF: G0 = 9.80665d0
A0 = 6375000.
G0 = 9.81D0
PI = 3.141592653589793D0
ONEPI = PI
PI180 = 180.D0/ONEPI
!---------GISS fix sigma-coords into hybrid eta-a and eta-b:
DO L=1,LTM+1
ETAA(L) = (1.E0 - SIGE(L))*PTOP
ETAB(L) = SIGE(L)
ENDDO
DO L=LTM+2,LM+1
ETAA(L) = SIGE(L)*PTOP
ETAB(L) = 0.0
ENDDO
!---------YGRD IS A SPECIAL CASE FOR REGULAR GISS GRID
YEDG0 = 180.D0 /REAL(JMG-1)
YDEDG(1) = -90. - YEDG0*0.5
DO J=2,JMG
YDEDG(J) = YDEDG(J-1) + YEDG0
ENDDO
YDEDG(1) = -90.
YDEDG(JMG+1) = +90
DO J=1,JMG+1
YEDG(J) = SIN(YDEDG(J)/PI180)
ENDDO
DO J=1,JMG
YGRD(J) = 0.5D0*(YEDG(J+1)+YEDG(J))
YDGRD(J) = 0.5D0*(YDEDG(J+1)+YDEDG(J))
ENDDO
DLONG = 2.D0*ONEPI/REAL(IMG)
XDEG0 = 0.D0
DO I=1,IMG
XGRD(I) = DLONG*REAL(I-1) + XDEG0/PI180
XDGRD(I) = PI180*XGRD(I)
IF(XDGRD(I).GT.180.D0) XDGRD(I) = XDGRD(I) - 360.D0
ENDDO
XEDG(1) = XGRD(1) - 0.5D0*DLONG
DO I=2,IMG
XEDG(I) = 0.5D0*(XGRD(I-1)+XGRD(I))
ENDDO
XEDG(IMG+1) = XGRD(IMG) + 0.5D0*DLONG
DO I=1,IMG+1
XDEDG(I) = PI180*XEDG(I)
IF(XDEDG(I).GT.180.D0) XDEDG(I) = XDEDG(I) - 360.D0
ENDDO
!---------areas:
DO J=1,JMG
DISTY(J) = A0*(YDEDG(J+1)-YDEDG(J))/PI180
ENDDO
DO J=1,JMG+1
DISTX(J) = A0*DLONG*SQRT(1.D0 - YEDG(J)*YEDG(J))
ENDDO
DO J=1,JMG
DISTX(J) = 0.5*(DISTX(J)+DISTX(J+1))
ENDDO
DO J=1,JMG
DO I=1,IMG
AREAXY(I,J) = A0*A0*(YEDG(J+1)-YEDG(J))*(XEDG(I+1)-XEDG(I))
ENDDO
ENDDO
!---------calculate air mass factors for each grid box:
!---------air mass (I,J,K) in kg = XYZA() + XYZB()*surface pressure (mbar)
SUMXYZ = 0.D0
DO K=1,KM
DO J=1,JM
JOFF = J+J0
DO I=1,IM
IOFF = MOD(I+I0-1,IMG)+1
!---set up for ECMWF inverse grid and Pascals (NOT mbar)
! XYZA(I,J,K) = (ETAA(K+1)-ETAA(K))*AREAXY(IOFF,JOFF)/G0
! XYZB(I,J,K) = (ETAB(K+1)-ETAB(K))*AREAXY(IOFF,JOFF)/G0
!---set for GISS grid, from bottom up, (will goto EC version)
XYZA(I,J,K) =-(ETAA(K+1)-ETAA(K))*AREAXY(IOFF,JOFF)*1.E2/G0
XYZB(I,J,K) =-(ETAB(K+1)-ETAB(K))*AREAXY(IOFF,JOFF)*1.E2/G0
SUMXYZ = SUMXYZ + XYZB(I,J,K)
ENDDO
ENDDO
ENDDO
DO J=1,JM
DO I=1,IM
XYA(I,J) = 0.D0
XYB(I,J) = 0.D0
DO K=1,KM
XYA(I,J) = XYA(I,J) + XYZA(I,J,K)
XYB(I,J) = XYB(I,J) + XYZB(I,J,K)
ENDDO
ENDDO
ENDDO
WRITE(lun13,'(A,1PE10.3)') ' Air Mass kg/mbar:',SUMXYZ
RETURN
END
!
!---------subroutines AIRSET-------------------------------------------
!
SUBROUTINE AIRSET 1
!---------AIRSET sets up the air mass-----------------------------------
REAL*8 AIRD,AIRNEW,AIRX, XYZA,XYZB,XYA,XYB
REAL*8 SUMXYZ,SUMAD0,SUMAW0,AIRWET,AIRH2O,AIRDRY,AIRQKG
!-----------------------------------------------------------------------
COMMON /CCGRID/ &
& AIRD(36,24,21),AIRNEW(36,24,21),AIRX(36,24,21), &
& XYZA(36,24,21),XYZB(36,24,21),XYA(36,24),XYB(36,24), &
& AIRDRY,AIRQKG, &
& YGRD(24),YDGRD(24),YEDG(25),YDEDG(25),XGRD(36), &
& XDGRD(36),XEDG(37),XDEDG(37),DLONG,AREAXY(36,24), &
& DISTX(25),DISTY(24), ETAA(22),ETAB(22), &
& UUU(37,24,21),VVV(36,25,21),PS(36,24), &
& AIRQG(36,24,21),AIRQ(36,24,21),PCTM(36,24), IMEPZ(24)
!-----------------------------------------------------------------------
# include "uci.h"
!xx$$
save
!xx$$
!-----------------------------------------------------------------------
!---------UNITS OF AIR MASS AND TRACER = (kg) -------------------------
!---------Air mass (kg) is given by area (m^2) * pressure thickness (Pa) / g0
!
! KWACC = KWACC+1
!
!---------GISS 21-layer, no q's, SET WATER MIXING RATIO SMALL EVERYWHERE:
IMG = IMX
JMG = JMX
KM = LM
G0 = 9.81D0
! write(6,*) 'Airset: airqg set to 0.d0'
DO J=1,JM
DO I=1,IM
DO K=1,KM
AIRQG(I,J,K) = 1.0D-6
! AIRQG(I,J,K) = 0.d0
ENDDO
ENDDO
ENDDO
!----------------------INITIALIZE---------------------------------------
!---------AIRD() = air mass at beginning of time step
!---------SUMAD0 = initialized dry-air mass (kg), this quantity
!--------- is assumed to be conserved throughout the run and is
!--------- used to reset AIRD()
SUMAD0 = 0.D0
SUMAW0 = 0.D0
!---------First need to map the full-grid values (__G) onto local(WINDOW) arrays
!---------NB need to have IMG=IM, JMG=JM, KMG=KM for full global calc.
DO J=1,JM
DO I=1,IM
II0 = MOD(I+I0-1, IMG) + 1
!---------original for q-code: PS(I,J) = PSG(II0,J+J0)--------------
PS(I,J) = P(II0,J+J0)
DO K=1,KM
AIRQ(I,J,K) = AIRQG(II0,J+J0,K)
ENDDO
ENDDO
ENDDO
!
DO K=1,KM
DO J=1,JM
DO I=1,IM
AIRWET = XYZA(I,J,K) + PS(I,J)*XYZB(I,J,K)
AIRH2O = AIRQ(I,J,K)*AIRWET
AIRD(I,J,K) = AIRWET - AIRH2O
SUMAD0 = SUMAD0 + AIRD(I,J,K)
SUMAW0 = SUMAW0 + AIRH2O
!---------temporary, duplication of AD() and AIRD() for current version
!rt AD(I,J,K) = AIRD(I,J,K)
!
ENDDO
ENDDO
ENDDO
!rt
AIRDRY = SUMAD0
IF (.NOT. LCONT) THEN
DO K=1,KM
DO J=1,JM
DO I=1,IM
AD(I,J,K) = AIRD(I,J,K)
ENDDO
ENDDO
ENDDO
ENDIF
!rt
!
! CALL TRASET
WRITE(lun13,'(A,1P,2E12.5)') ' Initialize Air: dry+water(kg):' &
& ,SUMAD0,SUMAW0
RETURN
END
!
!
SUBROUTINE PFILTR(XERR,AX,BX,NA,NB,NAM,NBM,NZONA,LSP,LNP,LEW) 1
!
IMPLICIT NONE
LOGICAL LSP,LNP,LEW
INTEGER NA,NB,NAM,NBM,NZONA(*)
REAL*8 XERR(NAM,NBM),AX(NAM+1,NBM),BX(NAM,NBM+1)
!-----------------------------------------------------------------------
! LSP = .T. if J=1 is S.POLE, LNP = .T. if J=NB is N.POLE
! LEW = .T. if I=1 connects to I=NA
! NZONA(J) = # of boxes in extended zones, offset by J0 when passed
!-----------------------------------------------------------------------
INTEGER I,IA,NAZ,J,J1,J2,NFLTR
REAL*8 X0(160,160)
REAL*8 SUMA,FNAZ8
!xx$$
save
!xx$$
!---------Initialize corrective column mass flows (kg): AX->alfa, BX->beta
DO J=1,NB+1
DO I=1,NA
BX(I,J) = 0.D0
ENDDO
ENDDO
DO J=1,NB
DO I=1,NA+1
AX(I,J) = 0.D0
ENDDO
ENDDO
DO J=1,NB
DO I=1,NA
X0(I,J) = XERR(I,J)
ENDDO
ENDDO
!---------Iterate over mass-error filter, accumulate corrections in AX & BX
DO 99 NFLTR=1,5
!---------must average XERR over EPZs, make sure that NZONA is offset by J0
DO J=1,NB
NAZ = NZONA(J)
IF(NAZ.GT.1) THEN
DO IA=1,NA,NAZ
SUMA = 0.D0
DO I=IA,IA+NAZ-1
SUMA = SUMA + XERR(I,J)
ENDDO
SUMA = SUMA/REAL(NAZ)
DO I=IA,IA+NAZ-1
XERR(I,J) = SUMA
ENDDO
ENDDO
ENDIF
ENDDO
!---------calculate BX = N-S filter, N-S wind between boxes (J-1) & (J)
DO J=3,NB-1
DO I=1,NA
BX(I,J) = BX(I,J) + 0.125D0*(XERR(I,J-1) - XERR(I,J))
ENDDO
ENDDO
!---------enhance the filtering by factor of 2 ONLY into/out-of polar caps
IF (LSP) THEN
DO I=1,NA
BX(I,2) = BX(I,2) + 0.25D0*(XERR(I,1) - XERR(I,2))
ENDDO
ELSE
DO I=1,NA
BX(I,1) = BX(I,1) - 0.125D0*XERR(I,1)
BX(I,2) = BX(I,2) + 0.125D0*(XERR(I,1) - XERR(I,2))
ENDDO
ENDIF
IF(LNP) THEN
DO I=1,NA
BX(I,NB) = BX(I,NB) + 0.25D0*(XERR(I,NB-1) - XERR(I,NB))
ENDDO
ELSE
DO I=1,NA
BX(I,NB+1) = BX(I,NB+1) + 0.125D0*XERR(I,NB)
BX(I,NB) = BX(I,J) + 0.125D0*(XERR(I,NB-1) - XERR(I,NB))
ENDDO
ENDIF
!---------need N-S flux across boundaries if window calc.(assume XERR=0 outside)
!---------NB for optimal matrix/looping, it would be best to define XERR=0
!---------for an oversized array XERR(0:IM+1,0:JM+1)
!---------calculate AX = E-W filter
J1 = 1
J2 = NB
IF (LSP) J1 = 2
IF (LNP) J2 = NB -1
DO J=J1,J2
!---------calculate pressure-filter E-W wind between boxes (I-1) & (I)
!---------enhance filtered wind by size of EPZ, will redistribute later within
NAZ = NZONA(J)
FNAZ8 = NAZ*0.125D0
DO I=2,NA
AX(I,J) = AX(I,J) + FNAZ8*(XERR(I-1,J) - XERR(I,J))
ENDDO
!---------calculate pressure-filter E-W wind at edges I=1 & I=NA+1
IF(LEW) THEN ! CYCLIC
AX(NA+1,J) = AX(NA+1,J) + FNAZ8*(XERR(NA,J) - XERR(1,J))
AX(1,J) = AX(1,J) + FNAZ8*(XERR(NA,J) - XERR(1,J))
ELSE ! WINDOW, assume zero error outside window
AX(1,J) = AX(1,J) - FNAZ8*XERR(1,J)
AX(NA+1,J) = AX(NA+1,J) + FNAZ8*XERR(NA,J)
ENDIF
!---------N.B. the AX()'s are NOT correct within EPZs, but ONLY at boundaries
!---------would need to fix below
! IF(NAZ.GT.1) THEN ! redistribute AX() at boundaries within EPZs
!
! ENDIF
ENDDO
!---------Update the mass error (XERR)
DO J=1,NB
DO I=1,NA
XERR(I,J) = X0(I,J) +AX(I,J) -AX(I+1,J) +BX(I,J) -BX(I,J+1)
ENDDO
ENDDO
99 CONTINUE
RETURN
END