C**** OIJV070.S make Atmospheric IxJ Vector datafile 1999/04/14 C**** C**** Compile: FCE90 OIJV070.S C**** C**** OIJV070 reads several Climate Model III accumulated diagnostic C**** files, calculates specified quantities from the OIJL arrays, C**** 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, * KOIJL0= JM*KAJ*6 + JM*LMA*KAJL + JM*3*4 + IM*JM*KAIJ + * IM*JM*LMA*KAIJL, * KDACC = KOIJL0 + 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$) C**** PARAMETER (KQMAX=3) LOGICAL*4 QQUAN(KQMAX), QMF,QEF,QSF INTEGER*4 IA(KQMAX), KOFCQ(2,KQMAX), LMINQ(KQMAX),LMAXQ(KQMAX), * LV(LMO) REAL*8 SCALE(KQMAX) CHARACTER FILEIN*80, TITLE*80, OUTYR*4,OUTMON*4, * NAME(KQMAX)*50, NAMEL*50, GRID(KQMAX)*1 COMMON /WORK00/ U(IM,JM),V(IM,JM) COMMON /QUANCB/ QMF,QEF,QSF, * LMINMF,LMINEF,LMINSF, LMAXMF,LMAXEF,LMAXSF EQUIVALENCE (QQUAN,QMF),(LMINQ,LMINMF),(LMAXQ,LMAXMF) C**** NAMELIST /INPUTZ/ LMINMF,LMINEF,LMINSF, LMAXMF,LMAXEF,LMAXSF, LV, * QMF, QEF, QSF DATA IA / 1 , 1 , 1 /, * KofCQ / 2,3, 6,7, 10,11/, * GRID / 'C', 'C', 'C'/ DATA SCALE /KQMAX*1./ DATA NAME / 1'OCEAN MASS FLUX (10^9 kg/s) of Layer 1', 2'POTENTIAL ENTHALPY FLUX (10^12 W) of Layer 1', 3'SALT MASS FLUX (10^6 kg/s) of Layer 1'/ DATA NAMEL /'OCEAN CURRENT (m/s) of Layer 1'/ C**** NARGS = IARGC() IF(NARGS.le.0) GO TO 800 C**** C**** Determine output quantities from Namelist input file C**** FILEIN = 'OIJV070.I ' OPEN (5,FILE=FILEIN,STATUS='OLD',ERR=810) READ (5,INPUTZ) CLOSE (5) OPEN (6,FILE='OIJV070.PRT') C**** C**** Zero out summing arrays C**** DO 10 N=1,15 10 IDACC(N) = 0 DO 20 I=1,IM*JM*LMO*KOIJL 20 OIJL(I,1,1,1) = 0. C**** C**** Loop over the input disk files C**** KFILE=0 100 KFILE=KFILE+1 IF(KFILE.gt.NARGS) GO TO 200 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,IHOURY CLOSE (1) IF(IHOURX.ne.IHOURY) GO TO 820 CALL STORYM (JYEAR0,JMON0, JYEAR,JMON) WRITE (6,911) JYEAR0,JMON0, JYEAR,JMON, FILEIN(1:60) C**** Accumulate diagnostics for time averaging DO 110 N=1,15 110 IDACC(N) = IDACC(N) + IDIAG(N) DO 120 I=1,IM*JM*LMO*KOIJL 120 OIJL(I,1,1,1) = OIJL(I,1,1,1) + DIAGR4(KOIJL0+I) GO TO 100 C**** C**** End of input files: print table of input years and months C**** 200 CALL PRNTYM (OUTYR,OUTMON) CALL SGEOM OPEN (2,FILE='OIJV070.O',FORM='UNFORMATTED') C**** Define scaling paramsters for some of the quantities SCALE(1) = 2.d- 9/NDYNO SCALE(2) = 1.d-12/DTS SCALE(3) = 1.d- 6/DTS C**** C**** Calculate vertically integrated diagnostics of mass flux (kg/s), C**** heat flux (W), and salt flux (kg/s) C**** DO 320 KQ=1,KQMAX IF(.not.QQUAN(KQ)) GO TO 320 U = 0. V = 0. KU = KofCQ(1,KQ) KV = KofCQ(2,KQ) DO 310 L=LMINQ(KQ),LMAXQ(KQ) DO 310 I=1,IM*JM U(I,1) = U(I,1) + OIJL(I,1,L,KU) 310 V(I,1) = V(I,1) + OIJL(I,1,L,KV) U = U*SCALE(KQ) / IDACC(IA(KQ)) V = V*SCALE(KQ) / IDACC(IA(KQ)) WRITE (NAME(KQ)(43:44),931) LMINQ(KQ) IF(LMAXQ(KQ).gt.LMINQ(KQ)) WRITE (NAME(KQ)(45:48),932) LMAXQ(KQ) CALL WRITED (NAME(KQ),LABEL,OUTMON,OUTYR,GRID(KQ)) 320 continue C**** C**** Ocean current (m/s) for individual layers C**** DO 430 K=1,LMO 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) = 4.*OIJL(I,J,L,2)*IDACC(5) / * ((OIJL(I,J,L,1)+OIJL(Ip1,J,L,1))*IDACC(1)*NDYNO*DYP(J)) 410 I=Ip1 DO 420 J=1,JM-1 DO 420 I=1,IM 420 V(I,J) = 4.*OIJL(I,J,L,3)*IDACC(5) / * ((OIJL(I,J,L,1)+OIJL(I,J+1,L,1))*IDACC(1)*NDYNO*DXV(J)) WRITE (NAMEL(29:30),931) L CALL WRITED (NAMEL,LABEL,OUTMON,OUTYR,'C') 430 continue C**** 500 CLOSE (2) CLOSE (6) GO TO 999 C**** 800 WRITE (0,*) 'Usage: OIJV070 /raid1/C070/DJan199* ', * 'make Ocean IxJ Vector files' WRITE (0,*) 'Check or edit input file: OIJV078.I 1999/04/14' GO TO 999 810 WRITE (0,*) 'Error accessing ',FILEIN STOP 810 820 WRITE (0,*) 'IHOURX and IHOURY do not match:',IHOURX,IHOURY STOP 820 C**** 911 FORMAT (' From',I6,A6,'to',I6,A6,' FILEIN=',A) 931 FORMAT (I2) 932 FORMAT (' -',I2) 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 OIJV070: inconsistant year and months received' WRITE (0,*) 'From OIJV070: 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,OUTMON,OUTYR,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,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,901) 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,902) TITLE, (U(IM/2+1,J),J=1,JM),(V(IM/2+1,J),J=1,JM) RETURN C**** 901 FORMAT (A50,A5,9X,A4,1X,A3) 902 FORMAT ('0',A72 / * 'U:',12F10.3 / 2X,12F10.3 / 2X,12F10.3 / 2X,10F10.3 / * 'V:',12F10.3 / 2X,12F10.3 / 2X,12F10.3 / 2X,10F10.3) END