C**** AIJV070.S make Atmospheric IxJ Vector datafile 98/04/08 C**** C**** AIJV070 reads several Climate Model III accumulated diagnostic C**** files, calculates specified quantities from the AIJ and AIJL C**** arrays, and writes the scaled data to a disk file. Each record C**** of the output file has the format: TITLE(80),U4(IM,JM),V4(IM,JM) C**** INCLUDE '/u/cmrun/C070.COM' PARAMETER (KDINT = 15 + 1 + 4 + 2*4 + 42*50, * KDACC = JM*KAJ*6 + JM*LMA*KAJL + JM*3*4 + IM*JM*KAIJ + * IM*JM*LMA*KAIJL + IM*JM*LMO*KOIJL + LMO*NMST*KOLNST + * 24*50*4 + JM*KCON) INTEGER*4 IPARM(300),IDIAG(KDINT) REAL*4 DIAGR4(KDACC) EQUIVALENCE (IPARM,IM$),(IDIAG,IDACC) C**** PARAMETER (KQMAX=8) LOGICAL*4 QQUAN(KQMAX), QSAV,QAMS,QSIMF,QSIV,QSIMS,QMF,QEF,QQF INTEGER*4 IDACCS(4), IA(KQMAX), KOFCQ(2,KQMAX), * LMINQ(6:8),LMAXQ(6:8), LV(LMA) REAL*8 U,V REAL*8 UVQ(IM,JM,KQMAX,2), SCALE(KQMAX) CHARACTER FILEIN*80, TITLE*80, NAME(KQMAX)*50, NAMEL*50, * OUTYR*4,OUTMON*4, GRID(KQMAX)*1 COMMON /WORK00/ U(IM,JM),V(IM,JM) COMMON /WORK01/ UQ(IM,JM,KQMAX),VQ(IM,JM,KQMAX),SRSI(IM,JM) COMMON /QUANCB/ QSAV,QAMS,QSIMF,QSIV,QSIMS, QMF,QEF,QQF, * LMINMF,LMINEF,LMINQF, LMAXMF,LMAXEF,LMAXQF EQUIVALENCE (UVQ,UQ),(QQUAN,QSAV), (LMINQ,LMINMF), (LMAXQ,LMAXMF) C**** NAMELIST /INPUTZ/ LMINMF,LMINEF,LMINQF,LMAXMF,LMAXEF,LMAXQF, LV, * QSAV, QAMS, QSIMF, QSIV, QSIMS, QMF, QEF, QQF DATA IA / 3 , 3 , 1 , 0 , 0 , 1 , 1 , 1 /, * KOFCQ /36,37, 48,49, 30,31, 67,68, 69,70, 2,3, 6,7, 10,11/, * GRID / 'A', 'A', 'C', 'C', 'C', 'C', 'C', 'C'/ DATA SCALE /KQMAX*1./ DATA NAME / 1'SURFACE AIR VELOCITY (m/s)', 2'ATMOSPHERIC MOMENTUM STRESS (N/m^2)', 3'SEA ICE MASS FLUX (10^6 kg/s)', 4'SEA ICE VELOCITY (m/s)', 5'SEA ICE MOMENTUM STRESS (N/m^2)', 6'ATMOSPHERIC MASS FLUX (10^9 kg/s) of Layer 1', 7'POTENTIAL ENTHALPY FLUX (10^15 W) of Layer 1', 8'WATER VAPOR MASS FLUX (10^6 kg/s) of Layer 1'/ DATA NAMEL /'ATMOSPHERIC VELOCITY (m/s) of Layer 1'/ C**** NARGS = IARGC() IF(NARGS.le.0) GO TO 800 OPEN (6,FILE='AIJV070.PRT') C**** Read first input file to determine KAIJ0, KAIJL0 and KDACCx CALL GETARG (1,FILEIN) OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) READ (1) IHOURX,LABEL,IPARM CLOSE (1) KAIJ0 = JM*KAJ$*6 + JM*LMA$*KAJL$ + JM*3*4 KAIJL0 = KAIJ0 + IM*JM*KAIJ$ KDACCx = KAIJL0 + IM*JM*LMA$*KAIJL$ + IM*JM*LMO$*KOIJL$ + + LMO$*NMST$*KOLNS$ + 24*50*4 + JM*KCON$ C**** Define scaling paramsters for some of the quantities SCALE(2) = NSURF/DTS SCALE(3) = 1.d-6/DTS SCALE(6) = 2.d-9/NDYNA SCALE(7) = 1.d-15*101325.**RKAP/DTS SCALE(8) = 1.d-6/DTS C**** C**** Determine output quantities from Namelist input file C**** FILEIN = 'AIJV070.I ' OPEN (5,FILE=FILEIN,STATUS='OLD',ERR=810) READ (5,INPUTZ) CLOSE (5) C**** Modify NAMEs DO 10 KQ=6,8 WRITE (NAME(KQ)(43:44),901) LMINQ(KQ) IF(LMAXQ(KQ).gt.LMINQ(KQ)) WRITE (NAME(KQ)(45:48),902) LMAXQ(KQ) 10 CONTINUE C**** C**** Zero out summing arrays C**** C UQ = VQ = WTQ = UCOM = VCOM = WTCOM = 0 C**** C**** Loop over the input disk files C**** DO 200 KFILE=1,NARGS CALL GETARG (KFILE,FILEIN) C**** Open and read input file OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) READ (1) IHOURX,LABEL,IPARM,IDIAG,(DIAGR4(K),K=1,KDACCx),IHOURY CLOSE (1) IF(IHOURX.ne.IHOURY) GO TO 811 CALL STORYM (JYEAR0,JMON0, JYEAR,JMON) WRITE (6,911) JYEAR0,JMON0, JYEAR,JMON, FILEIN(1:60) C**** C**** Sum accumulated diagnostics for time averaging, C**** DO 110 N=1,4 110 IDACCS(N) = IDACCS(N) + IDACC(N) DO 130 KQ=1,5 IF(.not.QQUAN(KQ)) GO TO 130 DO 120 KC=1,2 K = KOFCQ(KC,KQ) DO 120 I=1,IM*JM 120 UVQ(I,1,KQ,KC) = UVQ(I,1,KQ,KC) + DIAGR4(KAIJ0+IM*JM*(K-1)+I) 130 continue IF(QSIV.or.QSIMS) then DO 140 I=1,IM*JM 140 SRSI(I,1) = SRSI(I,1) + DIAGR4(KAIJ0+I) endif C**** Accumulate mass, energy and water vapor fluxes summed over layers DO 170 KQ=6,8 IF(.not.QQUAN(KQ)) GO TO 170 DO 160 KC=1,2 K = KOFCQ(KC,KQ) DO 160 L=LMINQ(KQ),LMAXQ(KQ) DO 160 I=1,IM*JM 160 UVQ(I,1,KQ,KC) = UVQ(I,1,KQ,KC) + + DIAGR4(KAIJL0+IM*JM*(L-1)+IM*JM*LMA$*(K-1)+I) 170 continue C**** Accumulate diagnostics for individual layer velocities DO 190 K=1,LMA IF(LV(K).le.0) GO TO 190 L =LV(K) DO 180 I=1,IM*JM MA(I,1,L) = MA(I,1,L) + DIAGR4(KAIJL0+IM*JM*(L-1)+I) UA(I,1,L) = UA(I,1,L) + DIAGR4(KAIJL0+IM*JM*(L-1)+IM*JM*LMA$+I) 180 VA(I,1,L) = VA(I,1,L) + DIAGR4(KAIJL0+IM*JM*(L-1)+IM*JM*LMA$*2+I) 190 continue 200 continue C**** C**** End of input files: print table of input years and months C**** CALL PRNTYM (OUTYR,OUTMON) OPEN (2,FILE='AIJV070.O',FORM='UNFORMATTED') C**** C**** Loop over quantities C**** DO 399 KQ=1,KQMAX IF(.not.QQUAN(KQ)) GO TO 399 GO TO (300,300,300,340,350,300,300,300),KQ C**** C**** Quantites that are scaled but not calculated C**** 300 DO 301 I=1,IM*JM U(I,1) = UQ(I,1,KQ)*SCALE(KQ) / IDACCS(IA(KQ)) 301 V(I,1) = VQ(I,1,KQ)*SCALE(KQ) / IDACCS(IA(KQ)) CALL WRITED (NAME(KQ),LABEL,OUTYR,OUTMON,GRID(KQ)) GO TO 399 C**** C**** Sea Ice Velocity (m/s) C**** 340 DO 341 I=1,IM SRSI(I,JM) = SRSI(1,JM) U(I,JM) = 0. 341 V(I,JM) = 0. I=IM DO 342 J=1,JM-1 DO 342 IP1=1,IM U(I,J) = 0. V(I,J) = 0. IF(SRSI(I,J)+SRSI(IP1,J).gt.0.) * U(I,J) = UQ(I,J,4) / (SRSI(I,J)+SRSI(IP1,J)) IF(SRSI(I,J)+SRSI(I,J+1).gt.0.) * V(I,J) = VQ(I,J,4) / (SRSI(I,J)+SRSI(I,J+1)) 342 I=IP1 CALL WRITED (NAME(4),LABEL,OUTYR,OUTMON,GRID(4)) GO TO 399 C**** C**** Sea Ice Momentum Stress (N/m^2) C**** 350 DO 351 I=1,IM SRSI(I,JM) = SRSI(1,JM) U(I,JM) = 0. 351 V(I,JM) = 0. I=IM DO 352 J=1,JM-1 DO 352 IP1=1,IM U(I,J) = 0. V(I,J) = 0. IF(SRSI(I,J)+SRSI(IP1,J).gt.0.) * U(I,J) = UQ(I,J,5) / (DTS*(SRSI(I,J)+SRSI(IP1,J))) IF(SRSI(I,J)+SRSI(I,J+1).gt.0.) * V(I,J) = VQ(I,J,5) / (DTS*(SRSI(I,J)+SRSI(I,J+1))) 352 I=IP1 CALL WRITED (NAME(5),LABEL,OUTYR,OUTMON,GRID(5)) 399 CONTINUE C**** C**** Atmospheric velocity (m/s) for individual layers C**** CALL SGEOM DO 430 K=1,LMA IF(LV(K).le.0) GO TO 430 L =LV(K) I=IM DO 410 J=1,JM DO 410 IP1=1,IM U(I,J) = 2.*IDACC(4)*UA(I,J,L) / / (NDYNA*IDACC(1)*.5*(MA(I,J,L)+MA(IP1,J,L))*DYP(J)) 410 I=IP1 DO 420 J=1,JM-1 DO 420 I=1,IM 420 V(I,J) = 2.*IDACC(4)*VA(I,J,L) / / (NDYNA*IDACC(1)*.5*(MA(I,J,L)+MA(IP1,J,L))*DYP(J)) WRITE (NAMEL(36:37),901) L CALL WRITED (NAMEL,LABEL,OUTYR,OUTMON,'C') 430 CONTINUE C**** 500 CLOSE (2) CLOSE (6) GO TO 999 C**** 800 WRITE (0,*) 'Usage: AIJV070 /raid1/C070/DJan198* ', * 'make Atmosphere IxJ Vector files' WRITE (0,*) 'Check or edit input file: AIJV070.I 98/04/08' GO TO 999 810 WRITE (0,*) 'Error accessing ',FILEIN STOP 810 811 WRITE (0,*) 'IHOURX and IHOURY do not match:',IHOURX,IHOURY STOP 811 C**** 901 FORMAT (I2) 902 FORMAT (' -',I2) 911 FORMAT (' From',I6,A5,' to',I6,A5,' FILEIN= ',A) 999 END SUBROUTINE SGEOM IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (IM=72,JM=46, * TWOPI=6.283185307179586477d0, RADIUS=6375000.) COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM) C**** C**** Calculates the spherical geometry for the C grid C**** DLON = TWOPI/IM DLAT = TWOPI*NINT(360./(JM-1))/720. FJEQ = .5*(1+JM) C**** Geometric parameters centered at secondary latitudes PLATS = -TWOPI/4. PSINS = SIN(PLATS) DO 10 J=1,JM-1 PLATN = DLAT*(J+1.-FJEQ) IF(J.eq.JM-1) PLATN = TWOPI/4. PSINN = SIN(PLATN) DXV(J) = DLON*RADIUS*(PSINN-PSINS) / (PLATN-PLATS) PLATS = PLATN 10 PSINS = PSINN C**** Geometric parameters centered at primary latitudes VLATS = -TWOPI/4. VSINS = SIN(VLATS) DO 20 J=1,JM VLATN = DLAT*(J+.5-FJEQ) IF(J.eq.JM) VLATN = TWOPI/4. VSINN = SIN(VLATN) DXYP(J)= DLON*RADIUS*RADIUS*(VSINN-VSINS) DXP(J) =.5*(DXV(J-1)+DXV(J)) DYP(J) = (VLATN-VLATS)*RADIUS VLATS = VLATN 20 VSINS = VSINN RETURN END SUBROUTINE STORYM (JYEAR0,JMON0, JYEAR,JMON) C**** C**** STORYM receives years and months and stores them in an array. C**** C**** Start of accumulation period: JYEAR0, JMON0, day 1, hour 0 C**** End of accumulation period: JYEAR, JMON, day 1, hour 0 C**** INTEGER*4 NYofM(12) CHARACTER*1 QMY(12,0:2500) CHARACTER*4 JMON0,JMON, MON(12), OUTYR,OUTMON DATA QMY /30012*' '/, NYofM/12*0/, MINYR/2501/, MAXYR/-1/ DATA MON /'Jan','Feb','Mar','Apr','May','Jun', * 'Jul','Aug','Sep','Oct','Nov','Dec'/ C**** IF(JMON0.eq.'Ann ') JMON0 = 'Jan ' ! compatible with prior vers IF(JYEAR0.gt.2500) JYEAR0 = 1900 + JYEAR0/100 DO 10 I=1,12 10 IF(JMON0.eq.MON(I)) GO TO 20 WRITE (6,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON0 is not a month.' WRITE (0,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON0 is not a month.' STOP 10 20 IMON0 = I DO 30 I=1,12 30 IF(JMON.eq.MON(I)) GO TO 40 WRITE (6,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON is not a month.' WRITE (0,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON is not a month.' STOP 30 40 IMON = I IF(IMON.gt.IMON0) then IMON = IMON-1 IYEAR = JYEAR else IMON = IMON+11 IYEAR = JYEAR-1 endif C**** DO 50 IY=JYEAR0,IYEAR DO 50 IM=IMON0 ,IMON IF(QMY(IM,IY).eq.'*') GO TO 805 QMY(IM,IY) = '*' JM = 1 + MOD(IM-1,12) NYofM(JM) = NYofM(JM) + 1 IF(JYEAR0.lt.MINYR) MINYR = JYEAR0 50 IF(JYEAR .gt.MAXYR) MAXYR = JYEAR RETURN C**** C**** C**** ENTRY PRNTYM (OUTYR,OUTMON) C**** C**** PRNTYM prints a table of input years and month, and C**** calculates OUTYR and OUTMON C**** C**** Output: OUTYR = character describing input years C**** OUTMON = character describing input months C**** C**** Produce printer table of months and years received C**** WRITE (6,910) DO 110 IY=MINYR,MAXYR 110 WRITE (6,911) IY,(QMY(IM,IY),IM=1,12) C**** C**** Determine output month: a single month, three months, or an annual C**** DO 210 M=1,12 210 IF(NYofM(M).gt.0) GO TO 220 GO TO 800 C**** Is output month a single month ? 220 M1 = M DO 230 M=M1+1,12 230 IF(NYofM(M).gt.0) GO TO 240 OUTMON = MON(M1) GO TO 400 C**** Is output twelve months ? 240 IF(M1.ne.1) GO TO 260 DO 250 M=2,12 250 IF(NYofM(M).ne.NYofM(M1)) GO TO 260 OUTMON = 'Ann ' GO TO 400 C**** Is output month three consecutive months ? 260 IF(M1.ge.11) GO TO 800 IF(NYofM(M1+1).ne.NYofM(M1)) GO TO 280 IF(NYofM(M1+2).ne.NYofM(M1)) GO TO 300 DO 270 M=M1+3,12 270 IF(NYofM(M).ne.0) GO TO 800 OUTMON = MON(M1)(1:1) // MON(M1+1)(1:1) // MON(M1+2)(1:1) // ' ' GO TO 400 C**** Is output 3 months November, December and January ? 280 IF(NYofM(11).ne.NYofM(1) .or. NYofM(12).ne.NYofM(1)) GO TO 800 DO 290 M=2,10 290 IF(NYofM(M).ne.0) GO TO 800 OUTMON = 'NDJ ' GO TO 400 C**** Is output 3 months December, January and February ? 300 IF(NYofM(12).ne.NYofM(1) .or. NYofM(2).ne.NYofM(1)) GO TO 800 DO 310 M=3,11 310 IF(NYofM(M).ne.0) GO TO 800 OUTMON = 'DJF ' C**** C**** Calculate the first and last years received, C**** and determine a character string describing those years C**** 400 MAXYR = MINYR + NYofM(M1) - 1 WRITE (OUTYR,940) MINYR IF(MINYR.lt.MAXYR) * WRITE (OUTYR,941) MOD(MINYR,100),MOD(MAXYR,100) IF(MINYR+8.le.MAXYR .and. MINYR/10.eq.MAXYR/10) * WRITE (OUTYR,942) MINYR/10 WRITE (6,943) OUTYR,OUTMON RETURN C**** 800 WRITE (6,*) 'From AIJV070: inconsistant year and months received' WRITE (0,*) 'From AIJV070: inconsistant year and months received' STOP 800 805 WRITE (6,901) JYEAR0,JMON0, JYEAR,JMON, 'accumulated already' WRITE (0,901) JYEAR0,JMON0, JYEAR,JMON, 'accumulated already' STOP 805 C**** 901 FORMAT (' From',I6,A5,' to',I6,A5,1X,A) 910 FORMAT ('0The following months and years were received:' / '0', * 9X,'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec') 911 FORMAT (1X,I6,21A5) 940 FORMAT (I4) 941 FORMAT (2I2.2) 942 FORMAT (I3,'X') 943 FORMAT ('0Description of the years and months above:',2A5) END SUBROUTINE WRITED (NAME,LABEL,OUTYR,OUTMON,GRID) C**** C**** Write a record to a data file whose format is: C**** TITLE(80),U(IM,JM),V(IM,JM) C**** C**** Input: U,V = vector components of quantity C**** NAME = name of quantity C**** LABEL = run number C**** OUTMON = month C**** OUTYR = year C**** GRID = location of components (for fixing poles) C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (IM=72,JM=46, TWOPI=6.283185307179586477d0) CHARACTER NAME*50,LABEL*6,OUTYR*4,OUTMON*4,GRID*1, TITLE*80 COMMON /WORK00/ U(IM,JM),V(IM,JM) C**** Rotate polar components for quantities on the A grid C**** (coding produces small error in V component at I=1) IF(GRID.ne.'A') GO TO 20 DO 10 I=IM,1,-1 ANGLE = TWOPI*(I-.5)/IM U(I, 1) = U(1, 1)*COS(ANGLE) - V(1, 1)*SIN(ANGLE) V(I, 1) = V(1, 1)*COS(ANGLE) + U(1, 1)*SIN(ANGLE) U(I,JM) = U(1,JM)*COS(ANGLE) + V(1,JM)*SIN(ANGLE) 10 V(I,JM) = V(1,JM)*COS(ANGLE) - U(1,JM)*SIN(ANGLE) C**** Build TITLE 20 WRITE (TITLE,902) NAME,LABEL,OUTYR,OUTMON C**** Write record to disk WRITE (2) TITLE,(SNGL(U(I,1)),I=1,IM*JM),(SNGL(V(I,1)),I=1,IM*JM) WRITE (6,903) TITLE,(U(IM/2+1,J),J=1,JM),(V(IM/2+1,J),J=1,JM) RETURN C**** 902 FORMAT (A50,A5,9X,A4,1X,A3) 903 FORMAT ('0',A72 / * ' U:',12F10.3 / 3X,12F10.3 / 3X,12F10.3 / 3X,10F10.3 / * ' V:',12F10.3 / 3X,12F10.3 / 3X,12F10.3 / 3X,10F10.3) END