C**** OTJ070.S North. Trans. of Mass, Heat and Salt 98/04/09 C**** C**** OTJ070 reads several Climate Model III D files, C**** calculates the northward transport of mass, heat and salt, C**** calculates the northward transport of heat by the horizontal gyre C**** and overturning, and writes the data as a Line Plot file. 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, * KOLNS0 = KOIJL0 + IM*JM*LMO*KOIJL, * KDACC = KOLNS0 + LMO*NMST*KOLNST + 24*50*4 + JM*KCON) INTEGER*4 IPARM(300),IDIAG(KDINT) REAL*4 DIAGR4(KDACC) EQUIVALENCE (IPARM,IM$),(IDIAG,IDACC) C**** CHARACTER FILEIN*80, OUTYR*4,OUTMON*4, * FILOUT(7)*20, TITLE*72, NAME(7)*50, YAXIS(3)*8 INTEGER*4 IS(4), KEYLOC(2,3) REAL*8 VLAT(0:JM),Q(0:JM,4,5),SCALE(3),SCALES(3),SOLNST(NMST), * MV(4),MT(4), XSPEC(3),YSPEC(3,3) COMMON /KBASCB/ KBASIN(IM,JM) C**** DATA FILOUT / 'OTJ070.MASS','OTJ070.HEAT','OTJ070.SALT', *'OTJ070.ATLN','OTJ070.PACF','OTJ070.INDN','OTJ070.GLOB'/ DATA NAME/ 1 'Northward Transport of Mass (10^9 kg/sec)', 2 'North. Trans. of Potential Enthalpy (10^15 W)', 3 'North. Trans. of Salt - .035*Mass (10^6 kg/sec)', 4 'North. Trans. of Heat in Atlantic Ocean (10^15 W)', 5 'North. Trans. of Heat in Pacific Ocean (10^15 W)', 6 'North. Trans. of Heat in Indian Ocean (10^15 W)', 7 'North. Trans. of Heat in Global Ocean (10^15 W)'/ DATA YAXIS /' Mass ', 'Enthalpy', ' Salt '/ DATA XSPEC /-90.,90.,30./ DATA YSPEC /-4.,4.,1., -5.,6.,1., -50.,50.,10./ DATA KEYLOC/401,350, 401,350, 401,350/ C**** C**** N Contents of OIJL(I,J,L,N) Scaling factor C**** - ------------------------- -------------- C**** 3 South-north mass flux (kg/s) 2/IDACC(1)*NDYNO C**** 7 South-north heat flux (J/s) 1/IDACC(1)*DTS C**** 11 South-north salt flux (kg/s) 1/IDACC(1)*DTS C**** C**** N Contents of OLNST(L,N) Scaling factor C**** - ---------------------- -------------- C**** 1 Strait mass flux (kg/s) 1/IDACC(1)*DTS C**** 2 Strait heat flux (J/s) .5/IDACC(1)*DTS C**** 3 Strait salt flux (kg/s) .5/IDACC(1)*DTS C**** NARGS = IARGC() IF(NARGS.le.0) GO TO 800 OPEN (6,FILE='OTJ070.PRT') C**** C**** Zero out summing arrays C**** IDACC1 = 0 DO 20 I=1,IM*JM*LMO*2 OIJL(I,1,1, 2) = 0. OIJL(I,1,1, 6) = 0. 20 OIJL(I,1,1,10) = 0. DO 30 L=1,LMO*NMST*3 30 OLNST(L,1,1) = 0. C**** C**** Loop over input D files, accumulate diagnostics C**** DO 150 KFILE=1,NARGS CALL GETARG (KFILE,FILEIN) C**** Open and read input D 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 811 CALL STORYM (JYEAR0,JMON0, JYEAR,JMON) WRITE (6,911) JYEAR0,JMON0, JYEAR,JMON, FILEIN(1:60) C**** Accumulate diagnostics for time averaging IDACC1 = IDACC1 + IDACC(1) DO 120 I=1,IM*JM*LMO*2 OIJL(I,1,1, 2) = OIJL(I,1,1, 2) + DIAGR4(KOIJL0+IM*JM*LMO +I) OIJL(I,1,1, 6) = OIJL(I,1,1, 6) + DIAGR4(KOIJL0+IM*JM*LMO*5+I) 120 OIJL(I,1,1,10) = OIJL(I,1,1,10) + DIAGR4(KOIJL0+IM*JM*LMO*9+I) DO 130 L=1,LMO*NMST*3 130 OLNST(L,1,1) = OLNST(L,1,1) + DIAGR4(KOLNS0+L) 150 continue C**** C**** End of input files: print table of input years and months C**** CALL PRNTYM (OUTYR,OUTMON) CALL OBASIN C**** NM = 3 NH = 7 SCALE(1) = 2.d- 9/(IDACC1*NDYNO) SCALE(2) = 1.d-15/(IDACC1*DTS) SCALE(3) = 1.d- 6/(IDACC1*DTS) SCALES(1) = 1.d- 9/(IDACC1*DTS) SCALES(2) = .5d-15/(IDACC1*DTS) SCALES(3) = .5d- 6/(IDACC1*DTS) C**** DO 190 KQ=1,5 DO 190 KB=1,4 DO 190 J=0,JM 190 Q(J,KB,KQ) = 0. C**** C**** Loop over quantites C**** DO 630 KQ=1,3 C**** C**** Calculate the transport entering each domain from the C**** southern boundary C**** NOIJL = 4*KQ-1 DO 220 J=1,JM-1 DO 220 I=1,IM IF(OIJL(I,J,1,NM).eq.0) GO TO 220 KB = KBASIN(I,J+1) SOIJL = 0. DO 210 L=1,LMO 210 SOIJL = SOIJL + OIJL(I,J,L,NOIJL) Q(J,KB,KQ) = Q(J,KB,KQ) + SOIJL*SCALE(KQ) Q(J, 4,KQ) = Q(J, 4,KQ) + SOIJL*SCALE(KQ) 220 continue C IF(KQ.ne.2) GO TO 300 C**** C**** Accumulate northward transport of heat by overturning and by C**** horizontal gyre C**** C DO 290 J=1,JM-1 C DO 280 L=1,LMO C DO 260 KB=1,3 C IS(KB) = 0 C MV(KB) = 0. C 260 MT(KB) = 0. C DO 270 I=1,IM C IF(OIJL(I,J,L,NM).eq.0) GO TO 270 C KB = KBASIN(I,J+1) C IS(KB) = IS(KB) + 1 C MV(KB) = MV(KB) + OIJL(I,J,L,NM) C MT(KB) = MT(KB) + OIJL(I,J,L,NH)/OIJL(I,J,L,NM) C 270 CONTINUE C IS(4) = IS(1) + IS(2) + IS(3) C MV(4) = MV(1) + MV(2) + MV(3) C MT(4) = MT(1) + MT(2) + MT(3) C DO 280 KB=1,4 C IF(IS(KB).eq.0) GO TO 280 C Q(J,KB,4) = Q(J,KB,4) + SCALE(2)*MV(KB)*MT(KB)/IS(KB) C 280 continue C DO 290 KB=1,4 C 290 Q(J,KB,5) = Q(J,KB,2) - Q(J,KB,4) C**** C**** Calculate the transport entering each domain from an east- C**** west boundary that is north of the southern boundary C**** C**** Bering Straight affects Atlantic and Pacific 300 SOIJL = 0. DO 310 L=1,LMO 310 SOIJL = SOIJL + OIJL(3,39,L,NOIJL) DO 320 J=1,39-1 Q(J,1,KQ) = Q(J,1,KQ) + SOIJL*SCALE(KQ) 320 Q(J,2,KQ) = Q(J,2,KQ) - SOIJL*SCALE(KQ) C**** Lombock Passage affects Indian and Pacific SOIJL = 0. DO 330 L=1,LMO DO 330 I=59,60 330 SOIJL = SOIJL + OIJL(I,21,L,NOIJL) DO 340 J=1,21-1 Q(J,2,KQ) = Q(J,2,KQ) + SOIJL*SCALE(KQ) 340 Q(J,3,KQ) = Q(J,3,KQ) - SOIJL*SCALE(KQ) C**** Borneo Passage affects Indian and Pacific SOIJL = 0. DO 350 L=1,LMO 350 SOIJL = SOIJL + OIJL(62,20,L,NOIJL) DO 360 J=1,20-1 Q(J,2,KQ) = Q(J,2,KQ) + SOIJL*SCALE(KQ) 360 Q(J,3,KQ) = Q(J,3,KQ) - SOIJL*SCALE(KQ) C**** C**** Calculate the transport entering each domain from a north- C**** south boundary that is north of the southern boundary C**** NOIJL = 4*KQ-2 C**** Drake Passage affects Pacific and Atlantic SOIJL = 0. DO 420 J=8,1,-1 DO 410 L=1,LMO 410 SOIJL = SOIJL + OIJL(23,J+1,L,NOIJL) Q(J,1,KQ) = Q(J,1,KQ) + SOIJL*SCALE(KQ) 420 Q(J,2,KQ) = Q(J,2,KQ) - SOIJL*SCALE(KQ) C**** Line south of Africa affects Atlantic and Indian SOIJL = 0. DO 440 J=13,1,-1 DO 430 L=1,LMO 430 SOIJL = SOIJL + OIJL(40,J+1,L,NOIJL) Q(J,1,KQ) = Q(J,1,KQ) - SOIJL*SCALE(KQ) 440 Q(J,3,KQ) = Q(J,3,KQ) + SOIJL*SCALE(KQ) C**** Line south of Austalia affects Indian and Pacific SOIJL = 0. DO 460 J=12,1,-1 DO 450 L=1,LMO 450 SOIJL = SOIJL + OIJL(65,J+1,L,NOIJL) Q(J,2,KQ) = Q(J,2,KQ) + SOIJL*SCALE(KQ) 460 Q(J,3,KQ) = Q(J,3,KQ) - SOIJL*SCALE(KQ) C**** C**** Calculate transport through straits from latitude to another C**** within the same basin C**** NOLNST=KQ DO 510 NS=1,NMST SOLNST(NS) = 0. DO 510 L=1,LMO 510 SOLNST(NS) = SOLNST(NS) + OLNST(L,NS,NOLNST) C**** Fury & Hecla: (19,42) to (20,40) Q(40,1,KQ) = Q(40,1,KQ) - SOLNST(1)*SCALES(KQ) Q(40,4,KQ) = Q(40,4,KQ) - SOLNST(1)*SCALES(KQ) Q(41,1,KQ) = Q(41,1,KQ) - SOLNST(1)*SCALES(KQ) Q(41,4,KQ) = Q(41,4,KQ) - SOLNST(1)*SCALES(KQ) C**** Nares: (22,43) to (24,44) Q(43,1,KQ) = Q(43,1,KQ) + SOLNST(2)*SCALES(KQ) Q(43,4,KQ) = Q(43,4,KQ) + SOLNST(2)*SCALES(KQ) C**** Gibralter: (35,32) to (37,33) Q(32,1,KQ) = Q(32,1,KQ) + SOLNST(3)*SCALES(KQ) Q(32,4,KQ) = Q(32,4,KQ) + SOLNST(3)*SCALES(KQ) C**** English: (36,36) to (37,37) Q(36,1,KQ) = Q(36,1,KQ) + SOLNST(4)*SCALES(KQ) Q(36,4,KQ) = Q(36,4,KQ) + SOLNST(4)*SCALES(KQ) C**** Bosporous: (42,33) to (43,34) Q(33,1,KQ) = Q(33,1,KQ) + SOLNST(6)*SCALES(KQ) Q(33,4,KQ) = Q(33,4,KQ) + SOLNST(6)*SCALES(KQ) C**** Red Sea: (44,29) to (45,28) Q(28,3,KQ) = Q(28,3,KQ) - SOLNST(7)*SCALES(KQ) Q(28,4,KQ) = Q(28,4,KQ) - SOLNST(7)*SCALES(KQ) C**** Bab-al-Mandab: (45,28) to (46,27) Q(27,3,KQ) = Q(27,3,KQ) - SOLNST(8)*SCALES(KQ) Q(27,4,KQ) = Q(27,4,KQ) - SOLNST(8)*SCALES(KQ) C**** Hormuz: (47,30) to (49,29) Q(29,3,KQ) = Q(29,3,KQ) - SOLNST(9)*SCALES(KQ) Q(29,4,KQ) = Q(29,4,KQ) - SOLNST(9)*SCALES(KQ) C**** Korea: (62,32) to (63,33) Q(32,2,KQ) = Q(32,2,KQ) + SOLNST(11)*SCALES(KQ) Q(32,4,KQ) = Q(32,4,KQ) + SOLNST(11)*SCALES(KQ) C**** Soya: (64,34) to (65,35) Q(34,2,KQ) = Q(34,2,KQ) + SOLNST(12)*SCALES(KQ) Q(34,4,KQ) = Q(34,4,KQ) + SOLNST(12)*SCALES(KQ) C**** C**** Calculate transport through straits from one basin to another C**** C**** Malacca: (56,25) to (58,24) DO 610 J=1,23 Q( J,2,KQ) = Q( J,2,KQ) + SOLNST(10)*SCALES(KQ) 610 Q( J,3,KQ) = Q( J,3,KQ) - SOLNST(10)*SCALES(KQ) Q(24,3,KQ) = Q(24,3,KQ) - SOLNST(10)*SCALES(KQ) Q(24,4,KQ) = Q(24,4,KQ) - SOLNST(10)*SCALES(KQ) C**** Fill in south polar value DO 620 KB=1,4 620 Q(0,KB,KQ) = Q(1,KB,KQ) 630 continue C**** Replace salt by salt - .035 times mass DO 640 KB=1,4 DO 640 J=0,JM 640 Q(J,KB,3) = Q(J,KB,3) - .035E3*Q(J,KB,1) C**** IDLAT= NINT(180./(JM-1)) FJEQ = JM/2. DO 660 J=1,JM-1 660 VLAT( J) = IDLAT*(J-FJEQ) VLAT( 0) = -90. VLAT(JM) = 90. C**** C**** Write titles and data to disk. C**** DO 710 KQ=1,3 OPEN (2,FILE=FILOUT(KQ)) WRITE (TITLE,970) NAME(KQ),LABEL(1:6),JMON0,JYEAR0 WRITE (6,*) ' ',TITLE(1:72) WRITE (2,*) TITLE WRITE (2,*) 'Latitude' WRITE (2,*) YAXIS(KQ) WRITE (2,*) ' Lat Atlantic Pacific Indian Global' WRITE (2,971) (VLAT(J),(Q(J,KB,KQ),KB=1,4),J=0,JM) C**** Information for plotting program WRITE (2,*) WRITE (2,972) ' X ', XSPEC WRITE (2,972) ' Y ',(YSPEC(K,KQ),K=1,3) WRITE (2,973) ' KEY',(KEYLOC(K,KQ),K=1,2) CLOSE (2) 710 continue C**** C**** Write titles and data to disk for northward transport of heat C**** by components C**** C DO 760 KB=1,4 C OPEN (2,FILE=FILOUT(3+KB)) C WRITE (TITLE,970) NAME(3+KB),LABEL(1:6),JMON0,JYEAR0 C WRITE (6,*) ' ',TITLE(1:72) C WRITE (2,*) TITLE C WRITE (2,*) 'Latitude' C WRITE (2,*) YAXIS(2) C WRITE (2,*) ' Lat Total Overturn Hor_Gyre' C WRITE (2,976) (VLAT(J),Q(J,KB,2),Q(J,KB,4),Q(J,KB,5),J=0,JM) C**** Information for plotting program C WRITE (2,*) C WRITE (2,972) ' X ', XSPEC C WRITE (2,972) ' Y ',(YSPEC(K,2),K=1,3) C WRITE (2,973) ' KEY',(KEYLOC(K,2),K=1,2) C CLOSE (2) C 760 continue GO TO 999 C**** 800 WRITE (0,*) 'Usage: OTJ070 /raid1/C070/D*1991 ', * 'North. transports by oceans 98/04/09' 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**** 911 FORMAT (' From',I6,A5,' to',I6,A5,' FILEIN=',A) 970 FORMAT (A50,' Run ',A6,2X,A3,I5) 971 FORMAT (F6.0,4F10.3) 972 FORMAT (A,1P,3E15.4) 973 FORMAT (A,2I6) 976 FORMAT (F6.0,3F10.3) 999 END SUBROUTINE OBASIN C**** C**** Read in KBASIN: 0=continent, 1=Atlantic, 2=Pacific, 3=Indian C**** PARAMETER (IM=72,JM=46) COMMON /KBASCB/ KBASIN(IM,JM) CHARACTER TITLE*72, CBASIN(IM,JM) C**** OPEN (3,FILE='/u/cmrun/KB4X510.OCN') READ (3,900) TITLE WRITE (6,*) 'Read on unit 3: ',TITLE READ (3,900) DO 10 J=JM,1,-1 10 READ (3,901) (CBASIN(I,J),I=1,72) CLOSE (3) DO 40 I=1,IM*JM KBASIN(I,1) = 0 IF(CBASIN(I,1).eq.'A') KBASIN(I,1) = 1 IF(CBASIN(I,1).eq.'P') KBASIN(I,1) = 2 IF(CBASIN(I,1).eq.'I') KBASIN(I,1) = 3 40 continue RETURN 900 FORMAT (A72) 901 FORMAT (72A1) 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 AIJK070: inconsistant year and months received' WRITE (0,*) 'From AIJK070: 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