1.   C****
   2.   C**** C088cM.S    Fortran source code for Main routines    2003/11/26
   3.   C****
   4.   C**** NASA/GISS Climate Model III               Gary L. Russell
   5.   C****                                           Reto A. Ruedy
   6.   C**** NASA/Goddard Space Flight Center          Gavin A. Schmidt
   7.   C**** Institute for Space Studies               Frank Abramopoulos
   8.   C**** 2880 Broadway, New York, NY 10025         Andrew A. Lacis
   9.   C**** U.S.A.                                    Jean A. Lerner
  10.   C****
  11.         PROGRAM MAIN
  12.         INCLUDE 'C070.COM'
  13.         PARAMETER (KDINT = 15 + 1 + 4 + 2*4 + 42*50,
  14.        *  KDACC = JM*KAJ*6 + JM*LMA*KAJL + JM*3*4 + IM*JM*KAIJ +
  15.        *          IM*JM*LMA*KAIJL + IM*JM*LMO*KOIJL + LMO*NMST*KOLNST +
  16.        *          24*50*4 + JM*KCON)
  17.         INTEGER*4 IPARM(300),IDIAG(KDINT)
  18.         REAL*8 ATMO(IM,JM,LMA,11),OCEN(IM,JM,LMO,11),GRND(IM,JM,4+LMG*4),
  19.        *       GICE(IM,JM,1+LMGI),SICE(IM,JM,8+LMSI),BERG(128,4),
  20.        *       RADN(IM,JM,9+LMA*2),STRA(LMO*7+4+LMSI,NMST), DIAG(KDACC)
  21.         EQUIVALENCE (IPARM,IM$),(ATMO,MA),(OCEN,MO),(GRND,SBG),
  22.        *  (GICE,MGI),(SICE,RSI),(BERG,BERGM),(STRA,MUST),(RADN,RQT),
  23.        *  (IDIAG,IDACC),(DIAG,AJ)
  24.   C****
  25.         PARAMETER (LMX=LMO, RHOI=910.)
  25.5        LOGICAL*4 QDIAGA
  26.         INTEGER*4 MSTART,MNOW,MSRCA
  30.         CHARACTER LABEND*4, MONYR*8, TITLE*80
  31.         COMMON /WORK00/ MM0(IM,JM,LMX),UM0(IM,JM,LMX),VM0(IM,JM,LMX)
  32.         COMMON /WORK01/ MM1(IM,JM,LMX),UM1(IM,JM,LMX),VM1(IM,JM,LMX)
  33.         COMMON /WORK02/ MMI(IM,JM,LMX),SMU(IM,JM,LMX),SMV(IM,JM,LMX),
  34.        *                SMW(IM,JM,LMO)
  36.         DATA SKIP /Z'FFEFFFFFFFFFFFFF'/
  37.   C****
  38.            CALL TIMER (MSTART,NDUMMY)
  39.   C**** Initialize common block parameters and prognostic variables
  40.         CALL INPUT
  41.            MSTART = MSTART-MDYNA-MDYNO-MCNDS-MRAD-MSURF-MSRCO-MDIAG-MELSE
  42.         CALL FFT0 (IM)
  45.              CALL CHECKT ('INPUT ')
  46.         QPK  = .FALSE.
  48.   C**** Initialize parameters used in MAIN
  49.         DTAFS = DTA*2./3.
  50.         DTALF = DTA*2.
  51.         DTOFS = DTO*2./3.
  52.         DTOLF = DTO*2.
  53.              CALL TIMER (MNOW,MELSE)
  54.              PERCNT = 100.*MELSE/(MNOW-MSTART)
  55.              WRITE (6,901) JYEAR,JMON,JDATE,JHOUR,MELSE,PERCNT,IHOUR
  56.         IF(IDT*24.ge.IHOURE*NDAY)  GO TO 630
  57.         WRITE (4,902) '####'
  58.         REWIND 4
  59.   C****
  60.   C**** Main loop
  61.   C****
  62.     100 IF(MOD(IDT-IDT0,NDISK).ne.0)  GO TO 120
  63.   C**** Write ReStartFile to disk
  64.         WRITE (KDISK) IHOUR,LABEL,IPARM,ATMO,OCEN,GRND,GICE,SICE,BERG,
  65.        *   STRA,RADN,BLDATA,IDIAG,DIAG,TMNMX,IHOUR
  66.         REWIND (KDISK)
  67.              CALL TIMER (MNOW,MELSE)
  68.              PERCNT = 100.*MELSE/(MNOW-MSTART)
  69.              WRITE (6,911) JYEAR,JMON,JDATE,JHOUR,MELSE,PERCNT,IHOUR,KDISK
  70.         KDISK = 3-KDISK
  71.   C**** Test for termination of run
  72.     120 READ (4,902) LABEND
  73.         REWIND (4)
  74.         IF(LABEND.eq.LABEL(1:4))  GO TO 800
  75.         IF(IDT*24.ge.IHOURE*NDAY)  GO TO 810
  76.   C****
  77.   C**** Initialize time parameters and diagnostic accumulating arrays
  78.   C**** if necessary (normally done at the beginning of a month)
  79.   C****
  80.            IF(IDT.le.IDT0)  GO TO 150
  81.            IF(MOD(IDT-IDT0,NDAY).ne.0)  GO TO 200
  82.            DO 140 K=0,12
  83.     140    IF(JDAY.eq.NDZERO(K))  GO TO 150
  84.            GO TO 180
  85.     150    JYEAR0 = JYEAR
  86.            JMON0  = JMON
  87.            JDATE0 = JDATE
  88.            JHOUR0 = JHOUR
  89.            JDAY0  = JDAY
  90.            DO 160 I=1,15
  91.     160    IDACC(I) = 0
  92.            IDACC(11) = 1
  93.            DO 170 K=1,KDACC
  94.     170    AJ(K,1) = 0.
  95.            CALL DIAGCA (1)
  96.            CALL DIAGCG (1)
  97.            CALL DIAGCL (1)
  98.            CALL DIAGCO (1)
  99.            CALL DIAGCI (1)
 100.            WRITE (6,918) JYEAR,JMON,JDATE,JHOUR,IHOUR
 101.   C**** Initialize some accumulating diagnostics at beginning of new day
 101.1    180    TMNMX(1,1,1) =  999.
 101.2           TMNMX(1,1,2) = -999.
 101.3           TMNMX(1,1,3) =  999.
 101.4           TMNMX(1,1,4) = -999.
 101.5           TMNMX(1,1,5) = 0.
 101.6           TMNMX(1,1,6) = 0.
 101.7           TMNMX(IM,1,7) = 0.
 101.8           DO 190 I=IM+1,IM*(JM-1)+1
 101.9           TMNMX(I,1,1) =  999.
 102.            TMNMX(I,1,2) = -999.
 102.1           TMNMX(I,1,3) =  999.
 102.2           TMNMX(I,1,4) = -999.
 102.3           TMNMX(I,1,5) = 0.
 102.4           TMNMX(I,1,6) = 0.
 102.5    190    TMNMX(I,1,7) = 0.
 106.            CALL TIMER (MNOW,MDIAG)
 109.   C****
 110.   C**** Integrate Atmospheric source terms
 111.   C****
 112.     200 QRAD = MOD(IDT-IDT0,NRAD).eq.0
 112.5             IDACC(1) = IDACC(1) + 1
 113.   C**** Condensation: moist convection and large scale
 114.         CALL MSTCNV
 115.         CALL CONDSE
 116.   C          CALL CHECKT ('CONDSE')
 118.         CALL PRECGI
 119.         CALL PRECOL
 120.         CALL PRECLI
 124.              CALL DIAGCA (4)
 124.1             CALL DIAGCG (4)
 124.2             CALL DIAGCL (4)
 125.5  C          CALL CHECKT ('PRECIP')
 125.6             CALL TIMER (MNOW,MCNDS)
 126.   C**** Radiation: solar and thermal
 127.         CALL RADIA
 128.   C          CALL CHECKT ('RADIA ')
 129.              CALL TIMER (MNOW,MRAD)
 130.              CALL DIAGCA (5)
 131.   C**** Surface interaction and ground calculation
 132.         CALL SURFCE
 133.   C          CALL CHECKT ('SURFCE')
 135.         CALL GLAICE
 136.         CALL OLAKE
 137.         CALL LAKICE
 139.   C          CALL CHECKT ('GROUND')
 140.         CALL DRYCNV
 142.              CALL DIAGCA (6)
 142.1             CALL DIAGCG (6)
 142.2             CALL DIAGCL (6)
 143.   C**** Stratospheric momentum drag
 144.         CALL ASDRAG
 144.1             CALL DIAGCA (7)
 144.5  C**** River Flow
 144.6        CALL RIVERF
 144.7             CALL DIAGCL (9)
 144.8             CALL DIAGCO (9)
 145.   C**** Replicate polar values to all longitudes
 146.         DO 210 L=1,LMA
 147.         DO 210 J=1,JM,JM-1
 148.         DO 210 I=2,IM
 149.         UA(I,J,L)  = UA(1,J,L)
 150.         H0M(I,J,L) = H0M(1,J,L)
 151.         HZM(I,J,L) = HZM(1,J,L)
 152.         Q0M(I,J,L) = Q0M(1,J,L)
 153.     210 QZM(I,J,L) = QZM(1,J,L)
 154.              CALL TIMER (MNOW,MSURF)
 156.   C          MSRCA  = MCNDS + MRAD + MSURF
 157.   C          PERCNT = 100.*MSRCA/(MNOW-MSTART)
 158.   C          WRITE (6,921) MSRCA,PERCNT
 159.   C****
 160.   C**** Integrate Ocean and Sea Ice source terms
 161.   C****
 162.         IF(.not.QOCEAN)  GO TO 300
 162.5        CALL KVINIT
 163.         CALL PRECOO
 164.         CALL PRECSI
 165.              CALL DIAGCO (4)
 166.              CALL DIAGCI (4)
 166.5        CALL DYNSI
 167.         CALL SEAICE
 167.5             CALL DIAGCO (5)
 167.6             CALL DIAGCI (5)
 168.         CALL OSTRES
 168.5             CALL DIAGCO (7)
 169.         CALL OSOURC
 169.5        CALL CREAIB
 170.              CALL DIAGCO (6)
 171.              CALL DIAGCI (6)
 172.   C          CALL CHECKT ('OSOURC')
 173.         CALL OCONV
 174.         CALL STCONV
 175.              CALL DIAGCO (8)
 176.         CALL OBDRAG
 177.         CALL OCOAST
 178.         CALL STDRAG
 179.              CALL DIAGCO (10)
 180.         CALL ADVSI
 181.   C     CALL STADVI
 182.         CALL MOVEIB
 183.              CALL DIAGCI (2)
 184.         CALL CALVEI
 185.              CALL DIAGCI (9)
 186.         CALL MELTSI
 187.         CALL MELTIB
 188.              CALL DIAGCO (11)
 189.              CALL DIAGCI (11)
 190.   C          CALL CHECKT ('MELTSI')
 193.   C**** Replicate polar values to all longitudes
 194.         DO 220 L=1,LMO
 195.                J=  JM
 196.         DO 220 I=2,IM
 197.         UO(I,J,L)  = UO(1,J,L)
 198.         G0M(I,J,L) = G0M(1,J,L)
 199.         GZM(I,J,L) = GZM(1,J,L)
 200.         S0M(I,J,L) = S0M(1,J,L)
 201.     220 SZM(I,J,L) = SZM(1,J,L)
 202.              CALL TIMER (MNOW,MSRCO)
 203.   C          PERCNT = 100*MSRCO/(MNOW-MSTART)
 204.   C          WRITE (6,922) MSRCO,PERCNT
 204.5        IF(.not.QATMO)  GO TO 400
 205.   C****
 206.   C**** Integrate Atmospheric Dynamic terms
 207.   C****
 209.     300 CALL AVTOM (MM0,UM0,VM0)
 210.         NS=NDYNA
 211.   C**** Initial forward step,  QMX = QM0 + .667*DTA*F(Q0)
 212.         CALL AFLUX  (NS,MM0,SKIP)
 213.         CALL AADVM (MM1,MM0,DTAFS)
 214.         CALL AADVV (UM1,VM1,UM0,VM0,DTAFS)
 215.         CALL APGF  (UM1,VM1,H0M,HZM,MM0,DTAFS)
 216.         CALL AMTOV (MM1,UM1,VM1)
 217.   C**** Initial backward step,  QM1 = QM0 + DTA*F(QX)
 218.         CALL AFLUX  (NS,MM0,SKIP)
 219.         CALL AADVM (MM1,MM0,DTA)
 220.         CALL AADVV (UM1,VM1,UM0,VM0,DTA)
 221.         CALL APGF  (UM1,VM1,H0M,HZM,MM0,DTA)
 222.         CALL AMTOV (MM1,UM1,VM1)
 223.         GO TO 340
 224.   C**** Odd leap frog step,  QM3 = QM1 + 2*DTA*F(Q2)
 225.     320 CALL AFLUX  (NS,MM1,SKIP)
 226.         CALL AADVM (MM1,MM1,DTALF)
 227.         CALL AADVV (UM1,VM1,UM1,VM1,DTALF)
 228.         CALL APGF  (UM1,VM1,H0M,HZM,MM0,DTALF)
 229.         CALL AMTOV (MM1,UM1,VM1)
 230.         NS=NS-1
 231.   C**** Even leap frog step,  QM2 = QM0 + 2*DTA*F(Q1)
 232.     340 CALL AFLUX  (NS,MM0,AIJL(1,1,1,2))
 233.         DO 350 I=1,IM*JM*LMA
 234.         MMI(I,1,1) = MM0(I,1,1)
 235.         SMU(I,1,1) = H0M(I,1,1)
 236.     350 SMV(I,1,1) = HZM(I,1,1)
 237.         CALL AADVT (H0M,HXM,HYM,HZM,MM0,DTALF,.FALSE.,AIJL(1,1,1,6))
 238.         CALL AADVT (Q0M,QXM,QYM,QZM,MM0,DTALF,.TRUE. ,AIJL(1,1,1,10))
 239.         CALL AADVM (MM0,MM0,DTALF)
 240.         CALL AADVV (UM0,VM0,UM0,VM0,DTALF)
 241.            CALL DIAGCD (1,DTALF)
 242.         DO 360 I=1,IM*JM*LMA
 243.         MMI(I,1,1) = .5*(MMI(I,1,1)+MM0(I,1,1))
 244.         SMU(I,1,1) = .5*(SMU(I,1,1)+H0M(I,1,1))  !  SMU = .5*(prior H0M +
 245.     360 SMV(I,1,1) = .5*(SMV(I,1,1)+HZM(I,1,1))  !            current H0M)
 246.         CALL APGF  (UM0,VM0,SMU,SMV,MMI,DTALF)
 247.            CALL DIAGCD (2,DTALF)
 247.5           QDIAGA = MOD((IDT-IDT0)*NDYNA-NS+2*(NDYNA/4),NDIAGA) .eq. 0
 248.            IF(QDIAGA)  CALL DIAGA
 249.         CALL AMTOV (MM0,UM0,VM0)
 250.   C**** Add prognostic variables to diagnostic AIJL arrays
 251.            IF(.not.QDIAGA)  GO TO 380
 252.            DO 370 I=1,IM*JM*LMA
 253.            AIJL(I,1,1,1) = AIJL(I,1,1,1) +  MA(I,1,1)
 254.            AIJL(I,1,1,5) = AIJL(I,1,1,5) + H0M(I,1,1)
 255.     370    AIJL(I,1,1,9) = AIJL(I,1,1,9) + Q0M(I,1,1)
 256.     380 NS=NS-1
 257.         IF(NS.gt.1)  GO TO 320
 258.   C****
 259.         QPK  = .FALSE.
 261.   C          CALL CHECKT ('ATMDYN')
 262.              CALL DIAGCA (2)
 263.              DO 390 I=IM,IM*(JM-1)+1
 264.              AIJ(I,1,14)  = AIJ(I,1,14)  + BLDATA(I,1,5)
 265.     390      TMNMX(I,1,7) = TMNMX(I,1,7) + BLDATA(I,1,5)
 266.   C**** Apply 8-th order Alternating Binomial Filter to UA and VA winds
 267.         CALL AABFIL
 267.1             CALL DIAGCA (3)
 267.5             CALL TIMER (MNOW,MDYNA)
 267.6  C          PERCNT = 100.*MDYNA/(MNOW-MSTART)
 267.7  C          WRITE (6,935) JYEAR,JMON,JDATE,JHOUR,MDYNA,PERCNT,IHOUR
 268.         IF(.not.QOCEAN)  GO TO 500
 269.   C****
 270.   C**** Integrate Ocean Dynamics terms
 271.   C****
 272.   C**** Add sea ice pressure to BLDATA(5) making BLDATA(6) (Pa-101325)
 273.     400 DO 405 I=IM+1,IM*(JM-1)+1
 274.     405 BLDATA(I,1,6) = BLDATA(I,1,5) +
 275.        +                GRAV*RSI(I,1)*(MSI(I,1,1)+MSI(I,1,2))
 277.         DO 410 I=1,IM*JM*(3*LMO-1)
 278.     410 SMU(I,1,1) = 0.
 279.         CALL OVTOM (MMI,UM0,VM0)
 280.         CALL OPGF0
 281.         NS=NDYNO
 282.   C**** Initial Forward step,  QMX = QM0 + DT*F(Q0)
 283.         CALL OFLUX  (NS,MMI,SKIP)
 284.         CALL OADVM (MM1,MMI,DTOFS)
 285.         CALL OADVV (UM1,VM1,UM0,VM0,DTOFS)
 286.         CALL OPGF  (UM1,VM1,DTOFS)
 287.         CALL OMTOV (MM1,UM1,VM1)
 288.   C**** Initial Backward step,  QM1 = QM0 + DT*F(Q0)
 289.         CALL OFLUX  (NS,MMI,SKIP)
 290.         CALL OADVM (MM1,MMI,DTO)
 291.         CALL OADVV (UM1,VM1,UM0,VM0,DTO)
 292.         CALL OPGF  (UM1,VM1,DTO)
 293.         CALL OMTOV (MM1,UM1,VM1)
 294.   C**** First even leap frog step,  Q2 = Q0 + 2*DT*F(Q1)
 295.         CALL OFLUX  (NS,MMI,SMU)
 296.         CALL OADVM (MM0,MMI,DTOLF)
 297.         CALL OADVV (UM0,VM0,UM0,VM0,DTOLF)
 298.         CALL OPGF  (UM0,VM0,DTOLF)
 299.         CALL OMTOV (MM0,UM0,VM0)
 300.         NS=NS-1
 301.   C**** Odd leap frog step,  Q3 = Q1 + 2*DT*F(Q2)
 302.     420 CALL OFLUX  (NS,MM1,SKIP)
 303.         CALL OADVM (MM1,MM1,DTOLF)
 304.         CALL OADVV (UM1,VM1,UM1,VM1,DTOLF)
 305.         CALL OPGF  (UM1,VM1,DTOLF)
 306.         CALL OMTOV (MM1,UM1,VM1)
 307.         NS=NS-1
 308.   C**** Even leap frog step,  Q4 = Q2 + 2*DT*F(Q3)
 309.         CALL OFLUX  (NS,MM0,SMU)
 310.         CALL OADVM (MM0,MM0,DTOLF)
 311.         CALL OADVV (UM0,VM0,UM0,VM0,DTOLF)
 312.         CALL OPGF  (UM0,VM0,DTOLF)
 313.         CALL OMTOV (MM0,UM0,VM0)
 314.         NS=NS-1
 315.         IF(NS.gt.1)  GO TO 420
 316.   C**** Advection of Potential Enthalpy and Salt
 317.            DO 430 I=1,IM*JM*(3*LMO-1)
 318.     430    OIJL(I,1,1,2) = OIJL(I,1,1,2) + SMU(I,1,1)
 319.         CALL OADVT (G0M,GXM,GYM,GZM,DTOLF,.FALSE.,OIJL(1,1,1,6))
 320.         CALL OADVT (S0M,SXM,SYM,SZM,DTOLF,.FALSE.,OIJL(1,1,1,10))
 321.              CALL CHECKT ('OCNDYN')
 321.5           IF(MOD(IDT-IDT0,NDIAGO).ne.0)  GO TO 450
 321.6           IDACC(5) = IDACC(5) + 1
 322.            DO 440 I=1,IM*JM*LMO
 323.            OIJL(I,1,1,1) = OIJL(I,1,1,1) +  MO(I,1,1)
 324.            OIJL(I,1,1,5) = OIJL(I,1,1,5) + G0M(I,1,1)
 325.     440    OIJL(I,1,1,9) = OIJL(I,1,1,9) + S0M(I,1,1)
 326.     450 QPK  = .FALSE.
 327.   C**** Advection of water through straits
 328.         CALL STPGF
 329.         CALL STADV
 330.              CALL DIAGCO (2)
 331.   C**** Apply Alternating Binomial Filter to UO and VO ocean currents
 332.         CALL OABFIL
 334.              CALL DIAGCO (3)
 334.1  C**** Add sea ice thickness to BLDATA(8) making BLDATA(7) (m^2/s^2)
 334.2        DO 470 I=IM+1,IM*(JM-1)+1
 334.3    470   AIJ(I,1,29) = BLDATA(I,1,8) + AIJ(I,1,29)
 334.4  C 470 BLDATA(I,1,7) = BLDATA(I,1,8) +
 334.5  C    +               (GRAV/RHOI)*RSI(I,1)*(MSI(I,1,1)+MSI(I,1,2))
 335.   C**** Replicate ocean and sea ice data to all longitudes for SIDYN
 336.         DO 480 I=2,IM
 337.   C     BLDATA(I,JM,5) = BLDATA(1,JM,5)
 338.         BLDATA(I,JM,6) = BLDATA(1,JM,6)
 338.1  C     BLDATA(I,JM,7) = BLDATA(1,JM,7)
 338.2    480 BLDATA(I,JM,8) = BLDATA(1,JM,8)
 339.              CALL TIMER (MNOW,MDYNO)
 340.   C          PERCNT = 100.*MDYNO/(MNOW-MSTART)
 341.   C          WRITE (6,945) JYEAR,JMON,JDATE,JHOUR,MDYNO,PERCNT,IHOUR
 342.   C****
 343.   C**** Update model time and call DAILY if required
 344.   C****
 345.     500 IDT   = IDT + 1
 346.         IHOUR = IDT*24/NDAY
 347.         IDAY  = IDT/NDAY
 348.         IYEAR = IDAY/365 + IYEAR0
 349.         JHOUR = MOD(IHOUR,24)
 350.         IF(MOD(IDT,NDAY).ne.0)  GO TO 510
 351.         CALL DAILY
 352.         CALL OCLIM  !  Climatology of lake temperatures and ice
 353.         CALL ACLIM  !  Atmospheric surface climatology
 353.5             CALL DIAGCA (12)
 353.6             CALL DIAGCL (12)
 354.              CALL TIMER (MNOW,MELSE)
 355.   C          PERCNT = 100.*MELSE/(MNOW-MSTART)
 356.   C          WRITE (6,945) JYEAR,JMON,JDATE,JHOUR,MELSE,PERCNT,IHOUR
 357.   C****
 358.   C**** Write special information to output files
 359.   C****
 381.     510 IF(NTAPE.eq.0.)  GO TO 600
 382.   C     IF(MOD(IDT-IDT0,NTAPE).ne.0)  GO TO 600
 383.   C     Computations for XXXXXX
 384.   C     WRITE (11) IHOUR,XXXXXX
 385.   C     ENDFILE 11
 386.   C     BACKSPACE 11
 387.   C          CALL TIMER (MNOW,MELSE)
 388.   C          PERCNT = 100.*MELSE/(MNOW-MSTART)
 389.   C          WRITE (6,959) MELSE,PERCNT,IHOUR
 390.   C****
 391.   C**** Call diagnostic printing routines
 392.   C****
 393.     600    IF(NDPRNT(0).ge.0)  GO TO 610
 394.            NDPRNT(0) = NDPRNT(0) + 1
 395.   C**** Print current diagnostics (used for the Initial Conditions)
 396.            IF(KDIAG(1).gt.0)  CALL DIAGCP
 397.            IF(KDIAG(2).gt.0)  CALL DIAGJ
 401.            CALL TIMER (MNOW,MDIAG)
 402.            GO TO 100
 403.   C**** Print diagnostic quantities on the NDPRNT-th day of the run
 404.     610    IF(MOD(IDT-IDT0,NDAY).ne.0)  GO TO 100
 405.            DO 620 K=0,12
 406.            IF(JDAY.eq.NDPRNT(K))  GO TO 630
 407.     620    continue
 408.            GO TO 700
 409.     630    IF(KDIAG(1).gt.0)  CALL DIAGCP
 410.            IF(KDIAG(2).gt.0)  CALL DIAGJ
 414.            IF(KDIAG(7).gt.0)  CALL DIAGD
 416.            CALL TIMER (MNOW,MDIAG)
 423.   C****
 424.   C**** Write output files at end of accumulation period
 425.   C****
 426.     700 DO 710 K=0,12
 426.1    710 IF(JDAY.eq.NDZERO(K))  GO TO 720
 426.2        GO TO 100
 427.   C**** Write model prognostic variables to output disk
 428.     720 WRITE (MONYR,970) 'M',JMON0,JYEAR0
 429.         OPEN  (3,FILE=MONYR,FORM='UNFORMATTED')
 430.         WRITE (3) IHOUR,LABEL,IPARM,ATMO,OCEN,GRND,GICE,SICE,BERG,STRA,
 431.        *  RADN,BLDATA,IHOUR
 432.         CLOSE (3)
 433.         WRITE (6,971) 'Model Variables',JYEAR,JMON,JDATE,JHOUR,IHOUR,MONYR
 434.   C**** Write accumulated diagnostics to output disk
 435.         WRITE (MONYR,970) 'D',JMON0,JYEAR0
 436.         OPEN  (3,FILE=MONYR,FORM='UNFORMATTED')
 437.         WRITE (3) IHOUR,LABEL,IPARM,IDIAG,(SNGL(DIAG(K)),K=1,KDACC),IHOUR
 438.         CLOSE (3)
 439.         WRITE (6,971) 'Acc Diagnostics',JYEAR,JMON,JDATE,JHOUR,IHOUR,MONYR
 440.         CALL TIMER (MNOW,MELSE)
 441.   C**** Print and zero out timing numbers at end of accumulation period
 442.         TOTALT= .01*(MNOW-MSTART)
 443.         PDYNA = MDYNA/TOTALT
 444.         PDYNO = MDYNO/TOTALT
 445.         PCDNS = MCNDS/TOTALT
 446.         PRAD  = MRAD /TOTALT
 447.         PSURF = MSURF/TOTALT
 448.         PSRCO = MSRCO/TOTALT
 449.         PDIAG = MDIAG/TOTALT
 450.         PELSE = MELSE/TOTALT
 451.         DTIME = (TOTALT/60.) /
 452.        *  ((JYEAR-JYEAR0)*365 + JDAY-JDAY0 + (JHOUR-JHOUR0)/24.d0)
 453.         WRITE (6,972) DTIME,PDYNA,PDYNO,PCDNS,PRAD,PSURF,PSRCO,PDIAG,PELSE
 454.         MDYNA = 0
 455.         MDYNO = 0
 456.         MCNDS = 0
 457.         MRAD  = 0
 458.         MSURF = 0
 459.         MSRCO = 0
 460.         MDIAG = 0
 461.         MELSE = 0
 462.         MSTART= MNOW
 464.         GO TO 100
 465.   C****
 466.   C**** End of main loop
 467.   C****
 468.   C**** Run terminated
 469.   C****
 470.     800 WRITE (6,980) LABEL(1:5)
 471.         WRITE (0,980) LABEL(1:5)
 472.     810 IF(MOD(IDT-IDT0,NDISK).eq.0)  GO TO 820
 473.         WRITE (KDISK) IHOUR,LABEL,IPARM,ATMO,OCEN,GRND,GICE,SICE,BERG,
 474.        *   STRA,RADN,BLDATA,IDIAG,DIAG,TMNMX,IHOUR
 475.         WRITE (6,981) JYEAR,JMON,JDATE,JHOUR,IHOUR,KDISK
 476.     820 WRITE (6,982) JYEAR,JMON,JDATE,JHOUR,IHOUR
 477.         WRITE (4,902) '####'
 478.         IF(IDT*24.ge.IHOURE*NDAY)  STOP 13
 479.         STOP 12
 480.   C****
 481.     901 FORMAT ('0NASA/GISS Climate Model III started     ',
 482.        *  I16,A5,I2,', Hr',I2,I16,F6.1,12X,'IHOUR=',I8)
 483.     902 FORMAT (A4)
 484.     911 FORMAT (' Restart file written to disk            ',
 485.        *  I16,A5,I2,', Hr',I2,I16,F6.1,12X,'IHOUR=',I8,'   on  fort.',I1/)
 486.     918 FORMAT (' Diagnostic accumulating arrays zeroed   ',
 487.        *  I16,A5,I2,', Hr',I2,34X,'IHOUR=',I8)
 488.     921 FORMAT (' Atmospheric source terms integrated     ',I45,F6.1)
 489.     922 FORMAT (' Ocean source terms integrated           ',I45,F6.1)
 490.     935 FORMAT (' Atmospheric dynamic terms integrated    ',
 491.        *  I16,A5,I2,', Hr',I2,I16,F6.1,12X,'IHOUR=',I8)
 492.     945 FORMAT (' Ocean dynamic terms integrated          ',
 493.        *  I16,A5,I2,', Hr',I2,I16,F6.1,12X,'IHOUR=',I8)
 494.     951 FORMAT (' Sea Level Pressure history written      ',
 495.        *  63X,'IHOUR=',I8,'   on  fort.10')
 496.     952 FORMAT (I3,'  U=',9F12.5)
 497.     959 FORMAT (' Special Information written             ',
 498.        *  I45,F6.1,12X,'IHOUR=',I8,'   on  fort.11')
 499.     970 FORMAT (A1,A3,I4)
 500.     971 FORMAT ('0',A15,' written to disk         ',
 501.        *  I16,A5,I2,', Hr',I2,34X,'IHOUR=',I8,'   on',A9/)
 502.     972 FORMAT (/'0Total time           Atmo Dyn   Ocean Dyn  ',
 503.        *     ' Condensat   Radiation     Surface   Ocean Src  ',
 504.        *     '  Diagnost       Other'/
 505.        * 1X,F5.2,' (min)    Percent:',F6.1,7F12.1)
 506.     980 FORMAT ('0Job  ',A5,'  terminated by END script.')
 507.     981 FORMAT (' ReStartFile written to disk             ',
 508.        *  I16,A5,I2,', Hr',I2,34X,'IHOUR=',I8,'   on  fort.',I1/)
 509.     982 FORMAT ('0NASA/GISS Climate Model III terminated normally',
 510.        *   I9,A5,I2,', Hr',I2,34X,'IHOUR=',I8)
 511.         END
 800.   
 801.         SUBROUTINE WRITE4 (TITLE,X1,X2,SCALE)
 802.   C****
 803.   C**** WRITE4 prints the contents of the REAL*4 array X
 804.   C****
 805.         PARAMETER (IM=72,JM=46, IMIN=31,JMIN=36,JMAX=40)
 806.         REAL*4 X1(IM,JM),X2(IM,JM),Y(10,JMAX-JMIN+1)
 807.         CHARACTER*16 TITLE
 808.         DO 10 J=JMIN,JMAX
 809.         DO 10 I=IMIN,IMIN+4
 810.         Y(I-IMIN+1,JMAX+1-J) = X1(I,J)/SCALE
 811.      10 Y(I-IMIN+6,JMAX+1-J) = X2(I,J)/SCALE
 812.         WRITE (6,901) TITLE,SCALE,Y
 813.         RETURN
 814.     901 FORMAT ('0',A16,'     units of',E15.6/(1X,5F13.3,2X,5F13.3))
 815.         END
 900.   
 901.         SUBROUTINE WRITE8 (TITLE,X1,X2,SCALE)
 902.   C****
 903.   C**** WRITE8 prints the contents of the REAL*8 array X
 904.   C****
 905.         PARAMETER (IM=72,JM=46, IMIN=31,JMIN=36,JMAX=40)
 906.         REAL*8 X1(IM,JM),X2(IM,JM),Y(10,JMAX-JMIN+1), SCALE
 907.         CHARACTER*16 TITLE
 908.         DO 10 J=JMIN,JMAX
 909.         DO 10 I=IMIN,IMIN+4
 910.         Y(I-IMIN+1,JMAX+1-J) = X1(I,J)/SCALE
 911.      10 Y(I-IMIN+6,JMAX+1-J) = X2(I,J)/SCALE
 912.         WRITE (6,901) TITLE,SCALE,Y
 913.         RETURN
 914.     901 FORMAT ('0',A16,'     units of',E15.6/(1X,5F13.3,2X,5F13.3))
 915.         END
1000.   
1001.         SUBROUTINE INPUT
1002.   C****
1003.   C**** INPUT reads in the starting atmospheric prognostic variables,
1004.   C**** overwrites some of the parameters in PARMCB, and
1005.   C**** reads in some of the fixed arrays
1006.   C****
1007.         INCLUDE 'C070.COM'
1008.         PARAMETER (KDINT = 16 + 4 + 2*4 + 42*50,
1009.        *  KDACC = JM*KAJ*6 + JM*LMA*KAJL + JM*3*4 + IM*JM*KAIJ +
1010.        *          IM*JM*LMA*KAIJL + IM*JM*LMO*KOIJL + LMO*NMST*KOLNST +
1011.        *          24*50*4 + JM*KCON)
1012.         INTEGER*4 IPARM(300),IDIAG(KDINT)
1013.         REAL*8 ATMO(IM,JM,LMA,11),OCEN(IM,JM,LMO,11),GRND(IM,JM,4+LMG*4),
1014.        *       GICE(IM,JM,1+LMGI),SICE(IM,JM,8+LMSI),BERG(128,4),
1015.        *       RADN(IM,JM,9+LMA*2),STRA(LMO*7+4+LMSI,NMST), DIAG(KDACC)
1016.         EQUIVALENCE (IPARM,IM$),(ATMO,MA),(OCEN,MO),(GRND,SBG),
1017.        *  (GICE,MGI),(SICE,RSI),(BERG,BERGM),(STRA,MUST),(RADN,RQT),
1018.        *  (IDIAG,IDACC),(DIAG,AJ)
1019.   C****
1020.         PARAMETER (RHOI=910., RHOS=100.)
1021.         REAL*4 VADATA
1022.         CHARACTER IDRUN*8, NLREC*80, TITLE*80, LABEL2*140
1023.         COMMON /RADNIO/ VADATA(11,4,3)
1024.         COMMON /WORK00/ SLP(IM,JM),Z500(IM,JM), JPARM(300)
1026.   C****
1027.         NAMELIST /INPUTZ/ ISTART, LS1,LSSM, QATMO,QOCEAN,QGRND,QCHECK,
1028.        *  IHOURI,IDAYI,IYEARI, IHOURE,IDAYE,IYEARE, IYEAR0,
1029.        *  NRAD,NSURF, NFILXA,NFILYA,NFILXO,NFILYO,
1030.        *  NDIAGO,NSLP,NTAPE,NDISK, IDEBUG,JDEBUG,
1031.        *  DTS,DTA,DTO, MDRYA,MTROP,MSTRAT,
1032.        *  SIGEA,SIGEO, IDACC,KEYCT,NAMDD,IJDD,KDIAG,NDZERO,NDPRNT
1033.         DATA ISTART /10/
1034.   C****
1035.         WRITE (6,901)
1036.         READ  (5,902) LABEL(1:72),LABEL(73:128)
1037.         WRITE (6,903) LABEL
1038.         IDRUN = LABEL(1:8)
1039.   C**** Copy INPUTZ namelist records into temporary disk file on unit 8
1040.         OPEN  (8,STATUS='SCRATCH')
1041.       5 READ  (5,904) NLREC
1042.         WRITE (8,904) NLREC
1043.         WRITE (6,905) NLREC
1044.         IF(NLREC(1:5).ne.' &END')  GO TO 5
1046.         IYEARI = IYEAR0
1047.         REWIND (8)
1048.         READ   (8,NML=INPUTZ)
1050.         IHOURI = ((IYEARI-IYEAR0)*365 + IDAYI)*24 + IHOURI
1051.   C****
1052.   C**** Calculates the spherical geometry for the C grid
1053.   C****
1054.         DLON = TWOPI/IM
1055.         DLAT = TWOPI*NINT(360./(JM-1))/720.
1056.   C**** Geometric parameters defined at secondary latitudes
1057.         FJEQ    = .5*(1+JM)
1058.         RLAT(1) = -TWOPI/4.
1059.         SINP(1) = SIN(RLAT(1))
1060.         DO 10 J=1,JM-1
1061.         RLAT(J+1)= DLAT*(J+1.-FJEQ)
1062.         IF(J.ge.JM-1)  RLAT(J+1) = TWOPI/4.
1063.         SINP(J+1)= SIN(RLAT(J+1))
1064.         COSV(J)  = (SINP(J+1)-SINP(J))/(RLAT(J+1)-RLAT(J))
1065.         DXV(J)   = DLON*RADIUS*COSV(J)
1066.         DYV(J)   = (RLAT(J+1)-RLAT(J))*RADIUS
1067.      10 DXYV(J)  = DLON*RADIUS*RADIUS*(SINP(J+1)-SINP(J))
1068.         DXV(JM)  = 0.
1069.         DYV(JM)  = 0.
1070.         DXYV(JM) = 0.
1071.         COSV(0)  = 0.
1072.         COSV(JM) = 0.
1073.   C**** Geometric parameters defined at primary latitudes
1074.         VLATS = -TWOPI/4.
1075.         SINVS = SIN(VLATS)
1076.         DO 20 J=1,JM
1077.         VLATN = DLAT*(J+.5-FJEQ)
1078.         IF(J.eq.JM)  VLATN = TWOPI/4.
1079.         SINVN    = SIN(VLATN)
1080.         DXYP(J)  = DLON*RADIUS*RADIUS*(SINVN-SINVS)
1081.         BYDXYP(J)= 1./DXYP(J)
1082.         DXP(J)   = .5*(DXV(J-1)+DXV(J))
1083.         DYP(J)   = (VLATN-VLATS)*RADIUS
1084.         COSP(J)  = .5*(COSV(J-1)+COSV(J))
1085.         DXYS(J)  = .5*DXYP(J)
1086.         DXYN(J)  = .5*DXYP(J)
1087.         VLATS    = VLATN
1088.      20 SINVS    = SINVN
1089.         DXP(1)   = .5*DXV(1)
1090.         DXP(JM)  = .5*DXV(JM-1)
1091.         DXYS(1)  = 0.
1092.         DXYS(JM) = DXYP(JM)
1093.         DXYN(1)  = DXYP(1)
1094.         DXYN(JM) = 0.
1095.   C**** Calculate area ratios for applying source terms to V winds
1096.         DO 30 J=1,JM-1
1097.         RAMVN(J  ) = DXYN(J  )/(DXYN(J)+DXYS(J+1))
1098.      30 RAMVS(J+1) = DXYS(J+1)/(DXYN(J)+DXYS(J+1))
1099.         RAMVS(1  ) = 0.
1100.         RAMVN(JM ) = 0.
1101.         DO 40 I=1,IM
1102.         SINI(I) = SIN((I-.5)*TWOPI/IM)
1103.      40 COSI(I) = COS((I-.5)*TWOPI/IM)
1104.   C****
1105.   C**** Set up vertical layering: calculate DSIGA and SIGA
1106.   C****
1107.         DO 50 L=1,LMA
1108.         DSIGA(L) = SIGEA(L-1)-SIGEA(L)
1109.         SIGA(L) = (SIGEA(L-1)+SIGEA(L))*.5
1110.      50 RDUMMY(L+49) = DSIGA(L)
1111.   C****
1112.   C**** Read in fixed arrays and calculate other arrays derived from them
1113.   C****
1114.   C**** Read in FIXDCB: FOCEAN, FLAKE, FGRND, FGICE, ZATMO and ZSOLID
1115.         CALL READR4 (11,IM*JM,FOCEAN,FOCEAN)
1116.         CALL READR4 (11,IM*JM,FLAKE ,FLAKE)
1117.         CALL READR4 (11,IM*JM,FGRND ,FGRND)
1118.         CALL READR4 (11,IM*JM,FGICE ,FGICE)
1119.         CALL READR4 (11,IM*JM,ZATMO ,ZATMO)
1120.         CALL READR4 (11,IM*JM,ZSOLID,ZSOLID)
1121.         CLOSE (11)
1122.   C**** Calculate FWATER and FLAND
1123.         DO 110 I=1,IM*JM
1124.         FWATER(I,1) = FOCEAN(I,1) + FLAKE(I,1)
1125.         FLAND(I,1)  = FGRND(I,1)  + FGICE(I,1)
1126.         IF(FWATER(I,1)+FLAND(I,1).ne.1.)  STOP 'INPUT 110'
1127.     110 continue
1128.   C**** Read in vegetation file: FVEG and VADATA
1129.         IF(.not.QGRND)  GO TO 150
1130.         DO 120 K=1,10
1131.     120 CALL READR4 (51,IM*JM,FVEG(1,1,K),FVEG(1,1,K))
1133.         READ  (51) TITLE,VADATA
1134.         CLOSE (51)
1135.   C**** Convert snow masking depth to (kg/m^2)
1136.         DO 130 K=1,10
1137.     130 VADATA(K,4,3) = VADATA(K,4,3)*RHOS
1138.         WRITE (6,912) (((VADATA(K,I,J),K=1,10),I=1,4),J=1,2),
1139.        *                (VADATA(K,4,3),K=1,10)
1140.   C**** Read in river direction file and fill in RIVRCB
1140.1        CALL RIVER0
1141.   C****
1142.     150 IF(ISTART.ge.10)  GO TO 600
1143.         GO TO (200,400,500),ISTART
1144.         GO TO 829
1145.   C****
1146.   C**** ISTART = 1:  Initialize run from initial conditions files
1147.   C****
1147.5    200 WRITE (6,920) ISTART,IHOURX,LABEL,
1147.6       *      'Model initialized from different I.C. files.'
1148.         IF(.not.QATMO)  GO TO 300
1148.1  C****
1148.2  C**** Read in and process atmospheric initial conditions from unit 10
1148.3  C****
1149.         CALL READR4 (10,IM*JM,MA,MA)
1149.1        DO 201 L=1,LMA*2
1149.2    201 CALL READR4 (10,IM*JM,UA(1,1,L),UA(1,1,L))
1149.3        DO 202 L=1,LMA
1149.4    202 CALL READR4 (10,IM*JM,H0M(1,1,L),H0M(1,1,L))
1149.5        DO 203 L=1,LMA
1149.6    203 CALL READR4 (10,IM*JM,Q0M(1,1,L),Q0M(1,1,L))
1150.         CLOSE (10)
1152.   C**** Set surface winds from first layer winds for radiation
1153.         BLDATA(1, 1,3) = SQRT(UA(1,2,1)*UA(1,2,1)+VA(1,1,1)*VA(1,1,1))
1154.         BLDATA(1,JM,3) = SQRT(UA(1,JM-1,1)*UA(1,JM-1,1)+
1155.        *                      VA(1,JM-1,1)*VA(1,JM-1,1))
1156.         IM1=IM
1157.         DO 210 J=2,JM-1
1158.         DO 210 I=1,IM
1159.         BLDATA(I,J,3) = .5*SQRT((UA(IM1,J,1)+UA(I,J,1))**2
1160.        *                       +(VA(I,J-1,1)+VA(I,J,1))**2)
1161.     210 IM1=I
1162.   C****
1163.         DO 270 J=1,JM
1164.         DO 270 I=1,IM
1165.   C**** Set radiation equilibrium heat from LMA layer temperature
1166.         DO 220 K=1,3
1167.     220 RQT(I,J,K) = H0M(I,J,LMA)*SHCD
1168.   C**** Replace temperature by potential temperature
1169.         DO 230 L=1,LMA
1170.     230 H0M(I,J,L) = H0M(I,J,L)/(SIGA(L)*100.*MA(I,J,1)+GRAV*MSTRAT)**RKAP
1170.5  C**** Store atmospheric surface pressure into BLDATA(5) (Pa-101325)
1170.6        BLDATA(I,J,5) = 100.*MA(I,J,1) + GRAV*MSTRAT - 101325.
1171.   C**** Convert Pressure (mb) to Mass per unit area (kg/m^2)
1172.         DO 240 L=LMA,1,-1
1173.     240 MA(I,J,L) = (100./GRAV)*DSIGA(L)*MA(I,J,1)
1175.   C**** Set stratospheric specific humidity to 3.d-6
1176.         DO 250 L=LS1,LMA
1177.     250 Q0M(I,J,L) = 3.d-6
1178.   C**** Define vertical potential temperature and humidity gradients
1179.         DO 260 L=2,LMA-1
1180.         HZM(I,J,L)  = .25*(H0M(I,J,L+1)-H0M(I,J,L-1))
1181.     260 QZM(I,J,L)  = .25*(Q0M(I,J,L+1)-Q0M(I,J,L-1))
1182.         HZM(I,J,1)  =  .5*(H0M(I,J,2)  -H0M(I,J,1))
1183.         QZM(I,J,1)  =  .5*(Q0M(I,J,2)  -Q0M(I,J,1))
1184.         HZM(I,J,LMA) = .5*(H0M(I,J,LMA)-H0M(I,J,LMA-1))
1185.         QZM(I,J,LMA) = .5*(Q0M(I,J,LMA)-Q0M(I,J,LMA-1))
1186.     270 continue
1187.   C**** Define horizontal gradients of potential temperature and humidity
1188.         DO 273 L=1,LMA
1189.         IM1=IM-1
1190.         I=IM
1191.         DO 271 J=1,JM-1
1192.         DO 271 IP1=1,IM
1193.         HXM(I,J,L) = .25*(H0M(IP1,J,L)-H0M(IM1,J,L))
1194.         HYM(I,J,L) = .25*(H0M(I,J+1,L)-H0M(I,J-1,L))
1195.         QXM(I,J,L) = .25*(Q0M(IP1,J,L)-Q0M(IM1,J,L))
1196.         QYM(I,J,L) = .25*(Q0M(I,J+1,L)-Q0M(I,J-1,L))
1197.         IM1=I
1198.     271 I=IP1
1199.         DO 272 J=1,JM,JM-1
1200.         DO 272 I=1,IM
1201.         HXM(I,J,L) = 0.
1202.         HYM(I,J,L) = 0.
1203.         QXM(I,J,L) = 0.
1204.     272 QYM(I,J,L) = 0.
1205.     273 continue
1206.   C**** Replace potential temperature by potential heat and
1207.   C**** specific humidity by water vapor mass
1208.         DO 280 L=1,LMA
1209.         DO 280 J=1,JM
1210.         DO 280 I=1,IM
1211.         H0M(I,J,L) = H0M(I,J,L)*(DXYP(J)*MA(I,J,L)*SHCD)
1212.         HXM(I,J,L) = HXM(I,J,L)*(DXYP(J)*MA(I,J,L)*SHCD)
1213.         HYM(I,J,L) = HYM(I,J,L)*(DXYP(J)*MA(I,J,L)*SHCD)
1214.         HZM(I,J,L) = HZM(I,J,L)*(DXYP(J)*MA(I,J,L)*SHCD)
1215.         Q0M(I,J,L) = Q0M(I,J,L)*(DXYP(J)*MA(I,J,L))
1216.         QXM(I,J,L) = QXM(I,J,L)*(DXYP(J)*MA(I,J,L))
1217.         QYM(I,J,L) = QYM(I,J,L)*(DXYP(J)*MA(I,J,L))
1218.     280 QZM(I,J,L) = QZM(I,J,L)*(DXYP(J)*MA(I,J,L))
1219.   C**** Define initial zonal velocities at the poles
1220.         DO 292 L=1,LMA
1221.         US = 0.
1222.         UN = 0.
1223.         DO 291 I=1,IM
1224.         US = US + UA(I,2,L)
1225.     291 UN = UN + UA(I,JM-1,L)
1226.         DO 292 I=1,IM
1227.         UA(I,1,L) = .25*US/IM
1228.     292 UA(I,JM,L)= .25*UN/IM
1229.   C****
1230.   C**** Read in ground initial conditions from unit 50
1231.   C**** Read in and process glacial ice initial conditions from unit 40
1232.   C****
1233.     300 IF(.not.QGRND)  GO TO 350
1234.         CALL READR4 (50,IM*JM,SBG,SBG)
1235.         DO 310 L=1,4*LMG+3
1236.     310 CALL READR4 (50,IM*JM,WBG(1,1,L),WBG(1,1,L))
1237.         CALL READR4 (50,IM*JM,BLDATA(1,1,1),BLDATA(1,1,1))
1238.         CALL READR4 (50,IM*JM,BLDATA(1,1,2),BLDATA(1,1,2))
1239.         CALL READR4 (50,IM*JM,BLDATA(1,1,4),BLDATA(1,1,4))
1240.         CLOSE (50)
1241.         CALL READR4 (40,IM*JM,MGI,MGI)
1242.         DO 320 L=1,LMGI
1243.     320 CALL READR4 (40,IM*JM,HGI(1,1,L),HGI(1,1,L))
1244.         CLOSE (40)
1245.   C**** Convert Glacial Ice temperatures into enthalpy minus latent heat
1246.         DO 330 I=1,IM*JM
1247.         HGI(I,1,1) = (SHCI*HGI(I,1,1)-ELHM)*XGI1*MGI(I,1)
1248.         HGI(I,1,2) = (SHCI*HGI(I,1,2)-ELHM)*XGI2*MGI(I,1)
1249.         HGI(I,1,3) = (SHCI*HGI(I,1,3)-ELHM)*XGI3*ACE2GI
1250.     330 HGI(I,1,4) = (SHCI*HGI(I,1,4)-ELHM)*XGI4*ACE2GI
1251.   C****
1252.   C**** Read in and process sea ice initial conditions from unit 30
1253.   C****
1254.     350 CALL READR4 (30,IM*JM,RSI,RSI)
1255.         CALL READR4 (30,IM*JM,MSI(1,1,1),MSI(1,1,1))
1256.         CALL READR4 (30,IM*JM,MSI(1,1,2),MSI(1,1,2))
1257.         DO 360 L=1,LMSI
1258.     360 CALL READR4 (30,IM*JM,HSI(1,1,L),HSI(1,1,L))
1259.         CLOSE (30)
1260.   C**** Zero out sea ice fraction over Ground or Glacial Ice
1261.         DO 370 I=1,IM*JM
1262.         IF(FLAND(I,1).ge.1.)  RSI(I,1) = 0.
1263.   C**** Convert RSI from a percentage to a ratio
1264.         RSI(I,1) = RSI(I,1)/100.
1268.   C**** Convert Sea Ice temperatures into enthalpy minus latent heat
1277.         HSI(I,1,1) = (SHCI*HSI(I,1,1)-ELHM)*XSI1*MSI(I,1,1)
1278.         HSI(I,1,2) = (SHCI*HSI(I,1,2)-ELHM)*XSI2*MSI(I,1,1)
1279.         HSI(I,1,3) = (SHCI*HSI(I,1,3)-ELHM)*XSI3*MSI(I,1,2)
1280.         HSI(I,1,4) = (SHCI*HSI(I,1,4)-ELHM)*XSI4*MSI(I,1,2)
1281.   C**** Load surface pressure below sea ice and surface geopotentials
1283.         BLDATA(I,1,6) = BLDATA(I,1,5) +
1284.        +                GRAV*RSI(I,1)*(MSI(I,1,1)+MSI(I,1,2))
1285.         BLDATA(I,1,7) = GRAV*ZATMO(I,1)
1286.         BLDATA(I,1,8) = BLDATA(I,1,7) -
1287.        -               (GRAV/RHOI)*RSI(I,1)*(MSI(I,1,1)+MSI(I,1,2))
1288.     370 continue
1289.         GO TO 450
1307.   C****
1308.   C**** ISTART = 2:  Initialize run from previous model output on unit 9.
1309.   C**** IPARM array is built from defaults and Namelist.
1310.   C****
1310.5    400 OPEN (9,FILE='fort.9',FORM='UNFORMATTED',STATUS='OLD')
1311.         READ (9,ERR=830) IHOURX,LABEL2,JPARM,ATMO,OCEN,GRND,GICE,SICE,
1312.        *     BERG,STRA,RADN,BLDATA
1313.         CLOSE (9)
1314.         WRITE (6,920) ISTART,IHOURX,LABEL2,
1315.        *      'Model initialized from previous model M file.'
1322.   C****
1323.   C**** New run is started: initialize model timing
1324.   C****
1325.     450 NDPRNT(0) = -1
1326.         IHOUR = IHOURI
1327.         IF(IHOURI.lt.0)  IHOUR = IHOURX
1327.5        NDAY  = NINT(SDAY/DTS)
1328.         IDT   = IHOUR*NDAY/24
1329.         IDAY  = IDT/NDAY
1330.         IYEAR = IDAY/365 + IYEAR0
1331.         JHOUR = MOD(IHOUR,24)
1332.         IDT0  = IDT
1333.         IHOUR0= IHOUR
1334.         IDAY0 = IDAY
1335.         CALL DAILY
1336.         CALL OCLIM
1337.         GO TO 650
1337.1  C****
1337.2  C**** ISTART = 3:  Restart Model from M file
1337.3  C****
1337.35   500 OPEN (9,FILE='fort.9',FORM='UNFORMATTED',STATUS='OLD')
1337.4        READ (9,END=830) IHOURX,LABEL,IPARM,ATMO,OCEN,GRND,GICE,SICE,
1337.5       *     BERG,STRA,RADN,BLDATA,IHOURY
1337.6        CLOSE (9)
1337.7        IF(IDRUN.ne.LABEL(1:8))  GO TO 864
1337.8        IF(IHOURX.ne.IHOURY)  GO TO 865
1337.9        WRITE (6,963) ISTART,IHOURX,LABEL,'M file.'
1338.         GO TO 650
1338.9  C****
1339.   C**** ISTART = 10-13:  Restart run from ReStartFiles 1 or 2
1340.   C****
1341.     600 KDISKX = ISTART-10
1342.         IF(ISTART.eq.11 .or. ISTART.eq.12)  GO TO 630
1343.   C**** Choose ReStartFile to restart from
1344.         IHOUR1 = -1
1345.         READ   (1,ERR=610) IHOUR1
1346.     610 REWIND (1)
1347.         IHOUR2 = -1
1348.         READ   (2,ERR=620) IHOUR2
1349.     620 REWIND (2)
1350.         IF(IHOUR1+IHOUR2.le.-2)  GO TO 862
1351.         KDISKX = 1
1352.         IF(IHOUR2.gt.IHOUR1)  KDISKX = 2
1353.         IF(ISTART.ge.13)  KDISKX = 3-KDISKX
1354.   C**** Restart on unit KDISKX
1355.     630 READ (KDISKX,ERR=863) IHOURX,LABEL,IPARM,ATMO,OCEN,GRND,GICE,
1356.        *  SICE,BERG,STRA,RADN,BLDATA,IDIAG,DIAG,TMNMX,IHOURY
1357.         REWIND (KDISKX)
1358.         IF(IDRUN.ne.LABEL(1:8))  GO TO 864
1359.         IF(IHOURX.ne.IHOURY)  GO TO 865
1360.         WRITE (6,963) ISTART,IHOURX,LABEL,'fort.',KDISK
1361.   C     KDISK = 3-KDISKX
1363.   C**** Update IPARM array from Namelist
1364.     650 IYEARE = IYEAR0
1365.         IDAYE  = 0
1366.         IHOURE = 0
1367.         REWIND (8)
1368.         READ   (8,NML=INPUTZ)
1369.         CLOSE  (8)
1370.         IHOURE = ((IYEARE-IYEAR0)*365 + IDAYE)*24 + IHOURE
1371.         LTM    = LS1-1
1372.         LMCM   = LTM
1372.5        NDAY   = NINT(SDAY/DTS)
1373.         NDYNA  = 2*NINT(.5*DTS/DTA)
1374.         NDYNO  = 2*NINT(.5*DTS/DTO)
1375.         DTA    = DTS/NDYNA
1376.         DTO    = DTS/NDYNO
1377.         NDIAGA = NDYNA*3+2
1378.         NDIAGS = NSURF+1
1378.5        CALL ACLIM
1379.   C****
1380.   C**** Reposition other files
1381.   C****
1382.   C**** Reposition output tape on unit 11
1383.   C     IF(NTAPE.le.0)  GO TO 720
1384.   C     KDISKX=11
1385.   C 710 READ (11,ERR=871,END=871)  IHOURX
1386.   C     IF(IHOUR*NDAY .ge. IHOURX*NDAY+NTAPE*24)  GO TO 710
1387.   C     WRITE (6,*) 'Output tape repositioned on unit 11.',
1388.   C    *  ' Last IHOURX read was: ',IHOURX
1389.   C 720 NTAPE=IABS(NTAPE)
1390.   C****
1391.         WRITE (6,INPUTZ)
1392.         CALL OINPUT (ISTART)
1393.         RETURN
1394.   C****
1395.   C**** Terminate because of improper start up
1396.   C****
1397.     829 WRITE (6,*) 'Invalid value for ISTART: ',ISTART
1398.         STOP 'INPUT 829'
1399.     830 WRITE (6,*) 'Error reading IC on unit 9.  IHOURX =',IHOURX
1400.         STOP 'INPUT 830'
1401.     831 WRITE (6,*) 'Error reading IC on unit 7.  IHOURY =',IHOURY
1402.         STOP 'INPUT 831'
1403.     862 STOP 'INPUT 862: Errors reading on both ReStartFiles.'
1404.     863 WRITE (6,*) 'Error reading ReStrartFile on unit:',KDISKX
1405.         IF(3-KDISKX.eq.KLAST)  GO TO 862
1406.         KLAST = KDISKX
1407.         KDISKX = 3-KDISK
1408.         WRITE (6,*) 'Attempt to read ReStartFile on unit:',KDISKX
1409.         GO TO 630
1410.     864 WRITE (6,*) 'Run Deck = ',IDRUN,'.  ReStartFile = ',LABEL(1:8)
1411.         STOP 'INPUT 864'
1412.     865 WRITE (6,*) 'IHOURX and IHOURY do not match on unit',KDISKX,
1413.        *  '.  IHOURX,IHOURY =',IHOURX,IHOURY
1414.         STOP 'INPUT 865'
1418.     871 WRITE (6,*) 'Error or EOF repositioning unit',KDISKX,
1419.        *  '.  Last IHOURX,IHOURY =',IHOURX,IHOURY
1420.         STOP 'INPUT 871'
1421.   C****
1422.     901 FORMAT (' NASA/GISS Climate Model III',30X,'1998 September 28',
1423.        *  20X,'Programmed by:  Gary L. Russell'/111X,'Reto A. Ruedy'/
1424.        *  ' NASA/Goddard Space Flight Center',78X,'Gavin A. Schmidt'/
1425.        *  ' Institute for Space Studies',83X,'Frank Abramopoulos'/
1426.        *  ' 2880 Broadway, New York, NY 10025',77X,'Andrew A. Lacis'/
1427.        *  ' U.S.A.',104X,'Jean A. Lerner'///)
1428.     902 FORMAT (A72/A56)
1429.     903 FORMAT ('0',A132/)
1430.     904 FORMAT (A80)
1431.     905 FORMAT (35X,A80)
1432.     912 FORMAT ('0Visual Albedoes:' / 10X,'Brite   Tndra   Grass   ',
1433.        *  'Shrub   Trees   Decid   Evrgr   RainF   Crops   Dark' /
1434.        *  ' Winter',10F8.4 / ' Spring',10F8.4 /
1435.        *  ' Summer',10F8.4 / ' Autumn',10F8.4 /
1436.        *       '0Near IR Albedoes:' / 10X,'Brite   Tndra   Grass   ',
1437.        *  'Shrub   Trees   Decid   Evrgr   RainF   Crops   Dark' /
1438.        *  ' Winter',10F8.4 / ' Spring',10F8.4 /
1439.        *  ' Summer',10F8.4 / ' Autumn',10F8.4 /
1440.        *'0Snow Masking Depth (kg/m^2):' / 10X,'Brite   Tndra   Grass   ',
1441.        *      'Shrub   Trees   Decid   Evrgr   RainF   Crops   Dark' /
1442.        *      7X,10F8.1)
1443.     920 FORMAT ('0ISTART =',I2,',  IHOUR =',I8 / 1X,A132 / 1X,A)
1444.     963 FORMAT ('0ISTART =',I2,',  IHOUR =',I8,',  Model restarted' /
1445.        *  3X,A80,'  from  ',A,I1)
1446.         END
1800.   
1801.         BLOCK DATA PARMBD
1802.   C****
1803.   C**** Default parameters for model common block
1804.   C****
1805.         INCLUDE 'C070.COM'
1806.         PARAMETER (p6=.6d0, p7=.7d0, p8=.8d0, mp6=-p6, mp7=-p7, mp8=-p8)
1810.         DATA IM$,JM$, LMA$,LMO$,LMSI$,LMGI$,LMG$, NMST$/
1811.        *     IM ,JM , LMA ,LMO ,LMSI ,LMGI ,LMG , NMST /,
1812.        *  KAJ$,KAJL$,KAIJ$,KAIJL$, KOIJL$,KOLNS$, KCON$/
1813.        *  KAJ ,KAJL ,KAIJ ,KAIJL , KOIJL ,KOLNST, KCON /,
1814.        *  NBERG, LS1,LSSM,  QATMO, QOCEAN, QGRND, QCHECK/
1815.        *      0,   8, LMA, .true.,.true. ,.true.,.false./,
1816.        *  NRAD,NSURF, NFILXA,NFILYA,NFILXO,NFILYO/
1817.        *     5,    2,      1,     0,     4,     6/,
1818.        *  NDIAGO,NSLP,NTAPE,NDISK, NDAY, KDISK, IYEAR0/
1819.        *       5,   0,    0,   24,   24,     1,   1958/,
1820.        *   DTS ,DTA ,DTO ,     MDRYA  ,   MTROP  ,  MSTRAT  /
1821.        *  3600.,450.,450., 10030.581d0,101.9368d0,101.9368d0/
1826.         DATA SIGEA /1.,.948665d0,.866530d0,.728953d0,.554415d0,
1827.        *               .390144d0,.251540d0,.143737d0,.061602d0,0./
1828.         DATA KEYCT /1/,
1829.        *  NAMDD  /'AUSD','MWST','SAHL','EPAC'/,
1830.        *  IJDD   / 63,17, 17,34, 37,26, 13,23/,
1831.        *  KEYNR  /2100*0/,
1832.        *  KDIAG  /1,1,14*0/,
1833.        *  NDZERO /0,1,32,60,91,121,152,182,213,244,274,305,335/,
1834.        *  NDPRNT /0,1,32,60,91,121,152,182,213,244,274,305,335/
1856.   C****
1857.   C****     Strait           From         To       LM    Width
1858.   C****     ------           ----         --       --    -----
1859.   C****  1  Fury & Hecla   19,42 ES    20,40 WN     2    20000
1860.   C****  2  Nares          22,43 EN    24,44 WS     5    30000 ! 50000
1861.   C****  3  Gibraltar      35,32 EN    37,33 WS     5    25000
1862.   C****  4  English        36,36 EN    37,37 WS     2    35000
1863.   C****  5  Kattegat       38,38 EN    40,38 WS     2    60000
1864.   C****  6  Bosporous      42,33 EN    43,34 WS     2     6000
1865.   C****  7  Red Sea        44,29 ES    45,28 WN     6   250000
1866.   C****  8  Bab al Mandab  45,28 ES    46,27 WN     6    25000
1867.   C****  9  Hormuz         47,30 ES    49,29 WN     2    50000
1868.   C**** 10  Malacca        56,25 EN    58,24 WS     3    50000
1869.   C**** 11  Korea          62,32 EN    63,33 WS     4   170000
1870.   C**** 12  Soya-kaikyo    64,34 EN    65,35 WS     2    40000
1871.   C****
1872.         DATA IST /19, 22, 35, 36, 38, 42, 44, 45, 47, 56, 62, 64,
1873.        *          20, 24, 37, 37, 40, 43, 45, 46, 49, 58, 63, 65/,
1874.        *     JST /42, 43, 32, 36, 38, 33, 29, 28, 30, 25, 32, 34,
1875.        *          40, 44, 33, 37, 38, 34, 28, 27, 29, 24, 33, 35/,
1876.        *     LMST/ 2,  5,  5,  2,  2,  2,  6,  6,  2,  3,  4,  2/,
1877.        *     XST / 0., p6, p6, 1., 1., 0., p6, p6, 1., p6, 0., p6,
1878.        *           0.,mp8,mp8, 0.,mp6,mp8,mp8,-1.,mp7,-1.,mp6,-1./,
1879.        *     YST /-1., p8, p8, 0., 0., 1.,mp8,mp8, 0.,mp8, 1., p8,
1880.        *           1.,mp6,mp6,-1.,mp8,mp6, p6, 0., p7, 0.,mp8, 0./,
1881.        *     WIST/ 20000., 30000., 25000., 35000., 60000.,  6000.,
1882.        *          250000., 25000., 50000., 50000.,170000., 40000./
1883.         END
1900.   
1901.         SUBROUTINE READR4 (IUNIT,NM,DATAR4,DATAR8)
1902.   C****
1903.   C**** READR4 reads a record from unit IUNIT containing TITLE,DATAR4,
1904.   C**** converts the REAL*4 array DATAR4 to the REAL*8 array DATAR8, and
1905.   C**** writes a line to unit 6 containing the TITLE just read.
1906.   C****
1907.         REAL*4 DATAR4(NM)
1908.         REAL*8 DATAR8(NM)
1909.         CHARACTER*80 TITLE
1910.   C****
1911.         READ (IUNIT) TITLE,DATAR4
1912.         DO 10 I=NM,1,-1
1913.      10 DATAR8(I) = DATAR4(I)
1914.         WRITE (6,901) IUNIT,TITLE
1915.         RETURN
1916.   C****
1917.     901 FORMAT (' Read on unit',I3,': ',A80)
1918.         END
2000.   
2001.         SUBROUTINE AVTOM (MM,UM,VM)
2002.   C****
2003.   C**** AVTOM converts density and velocities in concentration units
2004.   C**** to mass and momentum in mass units
2005.   C**** Input:  MA (kg/m**2), UA (m/s), VA (m/s)
2006.   C**** Output: MM (kg), UM (kg*m/s), VM (kg*m/s)
2007.   C****
2008.         IMPLICIT REAL*8 (A-H,M,O-Z)
2009.         PARAMETER (IM=72,JM=46,LMA=9)
2010.         REAL*8 MM(IM,JM,LMA),UM(IM,JM,LMA),VM(IM,JM,LMA)
2011.         COMMON /GEOMCB/ DXYP(JM)
2013.         COMMON /ATMOCB/ MA(IM,JM,LMA),UA(IM,JM,LMA),VA(IM,JM,LMA)
2014.         COMMON /WORK02/ MXY(IM,JM)
2014.5        DO 70 L=1,LMA
2015.   C**** Convert mass per unit area to mass
2016.         DO 10 J=1,JM
2017.         DO 10 I=1,IM
2018.      10 MM(I,J,L) = MA(I,J,L)*DXYP(J)
2019.   C**** Calculate mass weighting at B points
2020.         MXYS = 0.
2021.         MXYN = 0.
2022.         I=IM
2023.         DO 30 IP1=1,IM
2024.         DO 20 J=2,JM-2
2025.      20 MXY(I, J  ) = .125*(MM(I, J  ,L)+MM(IP1, J  ,L)+
2025.1       *                   (MM(I, J+1,L)+MM(IP1, J+1,L)))
2026.         MXY(I, 1  ) = .125*(MM(I, 2  ,L)+MM(IP1, 2  ,L)+MM(IM,1,L)*4.)
2027.         MXY(I,JM-1) = .125*(MM(I,JM-1,L)+MM(IP1,JM-1,L)+MM(1,JM,L)*4.)
2028.         MXYS = MXYS + MXY(I, 1  )
2029.         MXYN = MXYN + MXY(I,JM-1)
2030.      30 I=IP1
2031.   C**** Convert U velocity to U momentum
2033.         DO 40 J=2,JM-1
2034.         DO 40 I=1,IM
2035.      40 UM(I, J,L) = UA(I, J,L)*(MXY(I,J-1)+MXY(I,J))
2036.         UM(IM,1,L) = UA(IM,1,L)*MXYS
2037.         UM(1,JM,L) = UA(1,JM,L)*MXYN
2038.   C**** Convert V velocity to V momentum
2039.         IM1=IM
2040.         DO 60 J=1,JM-1
2041.         DO 60 I=1,IM
2042.         VM(I,J,L) = VA(I,J,L)*(MXY(IM1,J)+MXY(I,J))
2043.      60 IM1=I
2044.      70 continue
2045.         RETURN
2046.         END
2200.   
2201.         SUBROUTINE AMTOV (MM,UM,VM)
2202.   C****
2203.   C**** AMTOV converts mass and momentum in mass units
2204.   C**** to density and velocity in concentration units
2205.   C**** Input:  MM (kg), UM (kg*m/s), VM (kg*m/s)
2206.   C**** Output: MA (kg/m**2), UA (m/s), VA (m/s)
2207.   C****
2208.         IMPLICIT REAL*8 (A-H,M,O-Z)
2209.         PARAMETER (IM=72,JM=46,LMA=9)
2210.         REAL*8 MM(IM,JM,LMA),UM(IM,JM,LMA),VM(IM,JM,LMA)
2211.         COMMON /GEOMCB/ DXYP(JM)
2213.         COMMON /ATMOCB/ MA(IM,JM,LMA),UA(IM,JM,LMA),VA(IM,JM,LMA)
2214.         COMMON /WORK02/ MXY(IM,JM,LMA)
2214.5        DO 70 L=1,LMA
2215.   C**** Convert mass to mass per unit area
2216.         DO 10 J=1,JM
2217.         DO 10 I=1,IM
2218.      10 MA(I,J,L) = MM(I,J,L)/DXYP(J)
2219.   C**** Calculate mass weighting at B points
2220.         MXYS = 0.
2221.         MXYN = 0.
2222.         I=IM
2223.         DO 30 IP1=1,IM
2224.         DO 20 J=2,JM-2
2225.      20 MXY(I, J  ,L) = .125*(MM(I, J  ,L)+MM(IP1, J  ,L)+
2225.1       *                     (MM(I, J+1,L)+MM(IP1, J+1,L)))
2226.         MXY(I, 1  ,L) = .125*(MM(I, 2  ,L)+MM(IP1, 2  ,L)+MM(IM,1,L)*4.)
2227.         MXY(I,JM-1,L) = .125*(MM(I,JM-1,L)+MM(IP1,JM-1,L)+MM(1,JM,L)*4.)
2228.         MXYS = MXYS + MXY(I, 1  ,L)
2229.         MXYN = MXYN + MXY(I,JM-1,L)
2230.      30 I=IP1
2231.   C**** Convert U momentum to U velocity
2233.         DO 40 J=2,JM-1
2234.         DO 40 I=1,IM
2235.      40 UA(I, J,L) = UM(I, J,L)/(MXY(I,J-1,L)+MXY(I,J,L))
2236.         UA(1, 1,L) = UM(IM,1,L)/MXYS
2237.         UA(1,JM,L) = UM(1,JM,L)/MXYN
2238.         DO 50 I=2,IM
2239.         UA(I, 1,L) = UA(1, 1,L)
2240.      50 UA(I,JM,L) = UA(1,JM,L)
2241.   C**** Convert V momentum to V velocity
2242.         IM1=IM
2243.         DO 60 J=1,JM-1
2244.         DO 60 I=1,IM
2245.         VA(I,J,L) = VM(I,J,L)/(MXY(IM1,J,L)+MXY(I,J,L))
2246.      60 IM1=I
2247.      70 continue
2248.         RETURN
2249.         END
2400.   
2401.         SUBROUTINE AFLUX (NS,MM,AIJL)
2402.   C****
2403.   C**** AFLUX calculates the fluid fluxes
2404.   C**** Input:  MA (kg/m**2), UA (m/s), VA (m/s)
2405.   C**** Output: MU (kg/s) = DY * UA * MA
2406.   C****         MV (kg/s) = DX * VA * MA
2407.   C****         MW (kg/s) = MW(L-1) + CONV-SCONV*DSIGA + MM-SMM*DSIGA
2408.   C****       CONV (kg/s) = MU-MU+MV-MV
2409.   C****       AIJL (kg/s) = diagnostic accumulation of MU, MV and MW
2410.   C****
2411.         IMPLICIT REAL*8 (A-H,M,O-Z)
2412.         PARAMETER (IM=72,JM=46,LMA=9)
2413.         REAL*8 MM(IM,JM,LMA), AIJL(IM,JM,LMA,3)
2413.5        COMMON /PARMCB/ IPARM(100),DTS,DTA
2414.         COMMON /FLUXCB/ MU(IM,JM,LMA),MV(IM,JM,LMA),MW(IM,JM,LMA-1)
2415.         COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM),DYV(JM)
2416.         COMMON /LAYACB/ DSIGA(LMA)
2417.         COMMON /ATMOCB/ MA(IM,JM,LMA),UA(IM,JM,LMA),VA(IM,JM,LMA)
2418.         COMMON /WORK02/ DMVS(IM),DMVN(IM)
2419.         COMMON /WORK05/ CONV(IM,JM,LMA)
2420.         DATA SKIP /Z'FFEFFFFFFFFFFFFF'/
2421.   C****
2422.         DO 420 L=1,LMA
2423.   C**** Compute MU, the West-East mass flux, at non-polar points
2424.         I=IM
2425.         DO 120 J=2,JM-1
2426.         DO 120 IP1=1,IM
2427.         MU(I,J,L) = .5*DYP(J)*UA(I,J,L)*(MA(I,J,L)+MA(IP1,J,L))
2428.     120 I=IP1
2429.         CALL APFIL (MU(1,1,L))
2430.   C**** Compute MV, the South-North mass flux
2431.         DO 220 J=1,JM-1
2432.         DO 220 I=1,IM
2433.     220 MV(I,J,L) = .5*DXV(J)*VA(I,J,L)*(MA(I,J,L)+MA(I,J+1,L))
2434.   C**** Compute MU*2 at the poles
2435.         MUS = DYP( 1)*UA(1, 1,L)*MA(1, 1,L)
2436.         MUN = DYP(JM)*UA(1,JM,L)*MA(1,JM,L)
2437.         MVS = 0.
2438.         MVN = 0.
2439.         DO 310 I=1,IM
2440.         MVS = MVS + MV(I, 1  ,L)
2441.     310 MVN = MVN + MV(I,JM-1,L)
2442.         MVS = MVS/IM
2443.         MVN = MVN/IM
2444.         DMVS(1) = 0.
2445.         DMVN(1) = 0.
2446.         DO 320 I=2,IM
2447.         DMVS(I) = DMVS(I-1) + (MV(I, 1  ,L)-MVS)
2448.     320 DMVN(I) = DMVN(I-1) + (MV(I,JM-1,L)-MVN)
2449.         ADMVS = 0.
2450.         ADMVN = 0.
2451.         DO 330 I=1,IM
2452.         ADMVS = ADMVS + DMVS(I)
2453.     330 ADMVN = ADMVN + DMVN(I)
2454.         ADMVS = ADMVS/IM
2455.         ADMVN = ADMVN/IM
2456.         DO 340 I=1,IM
2457.         MU(I, 1,L) = (ADMVS-DMVS(I) + MUS)*2.
2458.     340 MU(I,JM,L) = (DMVN(I)-ADMVN + MUN)*2.
2468.   C****
2469.   C**** Compute horizontal fluid convergence
2470.   C****
2471.     400 IM1=IM
2472.         DO 410 J=2,JM-1
2473.         DO 410 I=1,IM
2474.         CONV(I,J,L) = MU(IM1,J,L)-MU(I,J,L) + (MV(I,J-1,L)-MV(I,J,L))
2475.     410 IM1=I
2476.         CONV(IM,1,L) = -MVS
2477.         CONV(1,JM,L) =  MVN
2478.     420 continue
2479.   C****
2480.   C**** Compute vertically integrated column convergence and mass
2481.   C****
2482.         DO 620 I=IM,IM*(JM-1)+1
2483.         SCONV = CONV(I,1,1)
2484.         SMM   =   MM(I,1,1)
2485.         DO 510 L=2,LMA
2485.1        SCONV = SCONV + CONV(I,1,L)
2485.2    510 SMM   =   SMM +   MM(I,1,L)
2486.   C****
2487.   C**** Compute MW, the upward fluid flux
2488.   C****
2489.         MW(I,1,1) = CONV(I,1,1)-SCONV*DSIGA(1) +
2489.1       *           (  MM(I,1,1)-  SMM*DSIGA(1))/(NS*DTA)
2490.         DO 610 L=2,LMA-1
2491.     610 MW(I,1,L) = CONV(I,1,L)-SCONV*DSIGA(L) + MW(I,1,L-1) +
2491.1       *           (  MM(I,1,L)-  SMM*DSIGA(L))/(NS*DTA)
2492.     620 continue
2493.   C****
2494.   C**** Accumulate diagnostics of MU, MV and MW
2495.   C****
2496.            IF(AIJL(1,1,1,1).eq.SKIP)  RETURN
2497.            DO 730 L=1,LMA
2498.            DO 710 I=IM+1,IM*(JM-1)
2499.     710    AIJL(I,1,L,1) = AIJL(I,1,L,1) + MU(I,1,L)
2500.            DO 720 I=1,IM
2501.            AIJL(I, 1,L,1) = AIJL(I, 1,L,1) + MU(I, 1,L)*.5
2502.     720    AIJL(I,JM,L,1) = AIJL(I,JM,L,1) + MU(I,JM,L)*.5
2503.            DO 730 I=1,IM*(JM-1)
2504.     730    AIJL(I,1,L,2) = AIJL(I,1,L,2) + MV(I,1,L)
2505.            DO 740 L=1,LMA-1
2506.            DO 740 I=IM,IM*(JM-1)+1
2507.     740    AIJL(I,1,L,3) = AIJL(I,1,L,3) + MW(I,1,L)
2508.         RETURN
2509.         END
2600.   
2601.         SUBROUTINE APFIL (X)
2602.   C****
2603.   C**** APFIL smoothes quantities in the zonal direction be reducing
2604.   C**** the coefficients of the Fourier series for high wave numbers
2605.   C**** near the poles.
2606.   C****
2606.5        IMPLICIT REAL*8 (A-H,M,O-Z)
2607.         PARAMETER (IM=72,JM=46,LMA=9, TWOPI=6.283185307179586477d0,
2607.5       *  DLON=TWOPI/IM, NMAX=IM/2)
2608.         INTEGER*4 NMIN(JM)
2609.         REAL*8 X(IM,JM),
2610.        *  BYSN(NMAX),DRAT(JM),SMOOTH(NMAX,JM),AN(0:NMAX),BN(0:NMAX)
2611.         COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM)
2612.   C****
2613.         DATA IFIRST /1/
2614.         IF(IFIRST.eq.0)  GO TO 100
2615.         IFIRST=0
2617.         DO 10 N=1,NMAX
2618.      10 BYSN(N) = 1./SIN(.5*DLON*N)
2619.         DO 30 J=2,JM-1
2620.         DRAT(J) = DXP(J)/DYP(3)
2621.         DO 20 N=NMAX,1,-1
2622.         SMOOTH(N,J) = BYSN(N)*DRAT(J)
2623.         IF(SMOOTH(N,J).GE.1.)  GO TO 30
2624.      20 continue
2625.   C     N=0
2626.      30 NMIN(J) = N+1
2627.   C****
2628.     100 DO 120 J=2,JM-1
2629.         IF(DRAT(J).GT.1)  GO TO 120
2630.         CALL FFT (X(1,J),AN,BN)
2631.         DO 110 N=NMIN(J),NMAX-1
2632.         AN(N) = AN(N)*SMOOTH(N,J)
2633.     110 BN(N) = BN(N)*SMOOTH(N,J)
2634.         AN(NMAX) = AN(NMAX)*SMOOTH(NMAX,J)
2635.         CALL FFTI (AN,BN,X(1,J))
2636.     120 continue
2637.         RETURN
2638.         END
2800.   
2801.         SUBROUTINE AADVM (MM2,MM0,DT)
2802.   C****
2803.   C**** AADVM calculates the updated column mass
2804.   C**** Input:  MM0 (kg), DT (s), CONV (kg/s), MW (kg/s)
2805.   C**** Output: MM2 (kg) = MM0 + DT*(CONV + MW(L-1)-MW(L))
2806.   C****
2807.         IMPLICIT REAL*8 (A-H,M,O-Z)
2808.         PARAMETER (IM=72,JM=46,LMA=9)
2809.         REAL*8 MM2(IM,JM,LMA),MM0(IM,JM,LMA)
2810.         COMMON /FLUXCB/ MU(IM,JM,LMA),MV(IM,JM,LMA),MW(IM,JM,LMA-1)
2811.         COMMON /GEOMCB/ DXYP(JM)
2812.         COMMON /ATMOCB/ MA(IM,JM,LMA),UA(IM,JM,LMA),VA(IM,JM,LMA)
2812.5        COMMON /LAYACB/ DSIGA(LMA)
2813.         COMMON /WORK05/ CONV(IM,JM,LMA)
2814.   C****
2815.   C**** Compute MM2, the new column mass
2816.   C****
2817.         DO 120 I=IM,IM*(JM-1)+1
2818.         MM2(I,1,1) = MM0(I,1,1) + DT*(CONV(I,1,1) - MW(I,1,1))
2818.1        DO 110 L=2,LMA-1
2818.2    110 MM2(I,1,L) = MM0(I,1,L) + DT*(CONV(I,1,L) + MW(I,1,L-1)-MW(I,1,L))
2818.3    120 MM2(I,1,LMA) = MM0(I,1,LMA) + DT*(CONV(I,1,LMA) + MW(I,1,LMA-1))
2819.   C**** Fill in values at the poles
2819.5        DO 130 L=1,LMA
2820.         DO 130 I=1,IM
2821.         MM2(I, 1,L) = MM2(IM,1,L)
2822.     130 MM2(I,JM,L) = MM2(1,JM,L)
2823.   C****
2824.   C**** Produce diagnostic output if new column mass exceeds limits
2825.   C****
2826.         DO 210 J=1,JM
2827.         IMAX=IM
2828.         IF(J.eq.1 .or. J.eq.JM)  IMAX=1
2829.         DO 210 I=1,IMAX
2830.         IF(MM2(I,J,1).GT.12000.*DSIGA(1)*DXYP(J))
2831.        *  WRITE (6,921) I,J,MM2(I,J,1)/DXYP(J),CONV(I,J,1),DT,
2832.        *    (UA(I-1,J,L),UA(I,J,L),VA(I,J-1,L),VA(I,J,L),L=1,LMA)
2833.     210 continue
2834.         RETURN
2835.   C****
2836.     921 FORMAT ('0Column mass exceeds limits.',2I5,F12.2,E15.5,F9.0/
2837.        *  ' UA(I-1,J,L)  UA(I,J,L)  VA(I,J-1,L)  VA(I,J,L)'/(4F11.4))
2838.         END
3000.   
3001.         SUBROUTINE AADVV (UM2,VM2,UM0,VM0,DT1)
3002.   C****
3003.   C**** AADVV advects atmospheric momentum (with coriolis force)
3004.   C**** Input:  MA (kg/m**2), UA (m/s), VA (m/s) = from odd solution
3005.   C****         MU (kg/s), MV (kg/s), MW (kg/s) = fluid fluxes
3006.   C**** Output: UM2 (kg*m/s) = UM0 + DT*(MU*U-MU*U+MV*U-MV*U+M*CM*V)
3007.   C****         VM2 (kg*m/s) = VM0 + DT*(MU*V-MU*V+MV*V-MV*V-M*CM*U)
3008.   C****
3008.5        IMPLICIT REAL*8 (A-H,M,O-Z)
3009.         PARAMETER (IM=72,JM=46,LMA=9, TWOPI=6.283185307179586477d0,
3010.        *  SDAY=86400., RADIUS=6375000., OMEGA=TWOPI*366./(365.*SDAY))
3011.         REAL*8 UM2(IM,JM,LMA),VM2(IM,JM,LMA),
3012.        *       UM0(IM,JM,LMA),VM0(IM,JM,LMA), DCOSP(JM),TANP(JM)
3013.         COMMON /FLUXCB/ MU(IM,JM,LMA),MV(IM,JM,LMA),MW(IM,JM,LMA-1)
3014.         COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM),DYV(JM),
3015.        *  COSV(0:JM),COSP(JM)
3018.         COMMON /ATMOCB/ MA(IM,JM,LMA),UA(IM,JM,LMA),VA(IM,JM,LMA)
3019.         COMMON /WORK03/ DUM(IM,JM,LMA),DVM(IM,JM,LMA),
3020.        *  FXY(IM,JM),GYY(IM,JM),GXY(IM,JM),FXX(IM,JM), UYGXYY(0:JM)
3021.         DATA IFIRST/1/
3022.   C****
3023.         IF(IFIRST.eq.0)  GO TO 20
3024.         IFIRST = 0
3025.         DO 10 J=1,JM
3026.         DCOSP(J) = COSV(J-1)-COSV(J)
3027.      10 TANP(J)  = DCOSP(J)/COSP(J)
3028.   C****
3029.      20 DT2  = DT1/2.
3030.         DT4  = DT1/4.
3031.         DT16 = DT1/16.
3032.         DT24 = DT1/24.
3033.         DT48 = DT1/48.
3034.   C**** Zero out momentum changes
3035.         DO 30 I=1,IM*JM*LMA*2
3036.      30 DUM(I,1,1) = 0.
3037.   C****
3038.   C**** Calculate combinations of MU and MV
3039.   C****
3040.         DO 480 L=1,LMA
3041.         IM1=IM-1
3042.         I=IM
3043.         DO 250 IP1=1,IM
3044.   C**** Calculate FXY, a 4 point average of MU, defined at V points
3045.         DO 210 J=1,JM-1
3046.     210 FXY(I,J) = MU(IM1,J,L)+MU(I,J,L)+(MU(IM1,J+1,L)+MU(I,J+1,L))
3047.   C**** Calculate GYY, a 3 point average of MV, defined at V points
3048.         DO 220 J=2,JM-2
3049.     220 GYY(I,J) = MV(I,J-1,L)+MV(I,J,L)+(MV(I,J,L)+MV(I,J+1,L))
3050.         GYY(I,1) = MV(I,1,L)+MV(I,2,L)
3051.         GYY(I,JM-1) = MV(I,JM-2,L)+MV(I,JM-1,L)
3052.   C**** Calculate GXY, a 4 point average of MV, defined at U points
3053.         DO 230 J=2,JM-1
3054.     230 GXY(I,J) = MV(I,J-1,L)+MV(IP1,J-1,L)+(MV(I,J,L)+MV(IP1,J,L))
3055.   C**** Calculate FXX, a 3 point average of MU, defined at U points
3056.         DO 240 J=2,JM-1
3057.     240 FXX(I,J) = MU(IM1,J,L)+MU(I,J,L)+(MU(I,J,L)+MU(IP1,J,L))
3058.         FXX(I,1) = 1.5*(MU(IM1, 1,L)+MU(I, 1,L)+(MU(I, 1,L)+MU(IP1, 1,L)))
3059.         FXX(I,JM)= 1.5*(MU(IM1,JM,L)+MU(I,JM,L)+(MU(I,JM,L)+MU(IP1,JM,L)))
3060.         IM1=I
3061.     250 I=IP1
3062.   C****
3063.   C**** Horizontal advection of momentum
3064.   C****
3065.   C     IM1=IM-1
3066.   C     I=IM
3067.         DO 350 IP1=1,IM
3068.   C**** Contribution of West-East flux to U wind
3069.         DO 310 J=2,JM-1
3070.         UFXYY = DT24*(FXY(I,J-1)+FXY(I,J))*(UA(IM1,J,L)+UA(I,J,L))
3071.         DUM(IM1,J,L) = DUM(IM1,J,L) - UFXYY
3072.     310 DUM(I  ,J,L) = DUM(I  ,J,L) + UFXYY
3073.   C**** Contribution of South-North and corner fluxes to U wind
3074.         DO 320 J=1,JM-1
3075.         UGXYY  = DT24*(GYY(I,J)+GYY(IP1,J))*(UA(I,J,L)+UA(I,J+1,L))
3076.         UFG45  = DT48*(GYY(I,J)+FXY(I,J))*(UA(IM1,J,L)+UA(I,J+1,L))
3077.         UFG135 = DT48*(GYY(I,J)-FXY(I,J))*(UA(I,J,L)+UA(IM1,J+1,L))
3078.         DUM(IM1,J  ,L) = DUM(IM1,J  ,L) - UFG45
3079.         DUM(I  ,J  ,L) = DUM(I  ,J  ,L) - (UGXYY+UFG135)
3080.         DUM(IM1,J+1,L) = DUM(IM1,J+1,L) + UFG135
3081.     320 DUM(I  ,J+1,L) = DUM(I  ,J+1,L) + (UGXYY+UFG45)
3082.   C**** Contribution of West-East flux to V wind
3083.         DO 330 J=1,JM-1
3084.         VFXXY = DT24*(FXX(I,J)+FXX(I,J+1))*(VA(I,J,L)+VA(IP1,J,L))
3085.         DVM(I  ,J,L) = DVM(I  ,J,L) - VFXXY
3086.     330 DVM(IP1,J,L) = DVM(IP1,J,L) + VFXXY
3087.   C**** Contribution of South-North and corner fluxes to V wind
3088.         DO 340 J=2,JM-1
3089.         VGXXY  = DT24*(GXY(IM1,J)+GXY(I,J))*(VA(I,J-1,L)+VA(I,J,L))
3090.         VFG45  = DT48*(GXY(I,J)+FXX(I,J))*(VA(I,J-1,L)+VA(IP1,J,L))
3091.         VFG135 = DT48*(GXY(I,J)-FXX(I,J))*(VA(IP1,J-1,L)+VA(I,J,L))
3092.         DVM(I  ,J-1,L) = DVM(I  ,J-1,L) - (VGXXY+VFG45)
3093.         DVM(IP1,J-1,L) = DVM(IP1,J-1,L) - VFG135
3094.         DVM(I  ,J  ,L) = DVM(I  ,J  ,L) + (VGXXY+VFG135)
3095.     340 DVM(IP1,J  ,L) = DVM(IP1,J  ,L) + VFG45
3096.         IM1=I
3097.     350 I=IP1
3098.   C****
3099.   C**** Coriolis force and metric term  (GXY is defined at the poles)
3100.   C****
3101.   C**** U component
3102.         UYGXYY(0) = 0.
3103.         UYGXYY(JM)= 0.
3104.         I=IM
3105.         DO 430 IP1=1,IM
3106.         GXY(I,1) = MV(I,1,L)+MV(IP1,1,L)
3107.         GXY(I,JM)= MV(I,JM-1,L)+MV(IP1,JM-1,L)
3108.         DO 410 J=1,JM-1
3109.     410 UYGXYY(J) = .0625*(UA(I,J,L)+UA(I,J+1,L))*(GYY(I,J)+GYY(IP1,J))
3110.         DO 420 J=1,JM
3111.     420 DUM(I,J,L) = DUM(I,J,L) + DT2*(OMEGA*RADIUS*GXY(I,J)*DCOSP(J) +
3112.        *  (UYGXYY(J-1)+UYGXYY(J))*TANP(J))
3113.     430 I=IP1
3114.   C**** V component
3115.         DO 450 I=1,IM
3116.         DO 440 J=1,JM-1
3117.     440 UYGXYY(J) = .0625*(UA(I,J,L)+UA(I,J+1,L))*
3118.        *  (UA(I,J,L)*TANP(J)+UA(I,J+1,L)*TANP(J+1))
3119.         FXY(I,1) = OMEGA*RADIUS*UA(1, 1,L)*DCOSP( 1)
3120.         FXY(I,JM)= OMEGA*RADIUS*UA(1,JM,L)*DCOSP(JM)
3121.         DO 450 J=2,JM-1
3122.     450 FXY(I,J) = OMEGA*RADIUS*UA(I,J,L)*DCOSP(J)+(UYGXYY(J-1)+UYGXYY(J))
3123.         IM1=IM
3124.         DO 470 I=1,IM
3125.         DO 460 J=1,JM-1
3126.     460 DVM(I,J,L) = DVM(I,J,L) - DT4*DXV(J)*(MA(I,J,L)+MA(I,J+1,L))*
3127.        *  (FXY(IM1,J)+FXY(I,J)+(FXY(IM1,J+1)+FXY(I,J+1)))
3128.     470 IM1=I
3129.     480 continue
3130.   C****
3131.   C**** Vertical advection of momentum
3132.   C****
3133.         DO 560 L=1,LMA-1
3134.         FXYS = 0.
3135.         FXYN = 0.
3136.         I=IM
3137.         DO 520 IP1=1,IM
3138.         DO 510 J=2,JM-2
3139.     510 FXY(I,J) = (MW(I,J,L)+MW(IP1,J,L)+(MW(I,J+1,L)+MW(IP1,J+1,L)))
3140.         FXY(I,1) = (MW(I,2,L)+MW(IP1,2,L)+4.*MW(IM,1,L))
3141.         FXY(I,JM-1)= (MW(I,JM-1,L)+MW(IP1,JM-1,L)+4.*MW(1,JM,L))
3142.         FXYS = FXYS + FXY(I,1)
3143.         FXYN = FXYN + FXY(I,JM-1)
3144.     520 I=IP1
3145.   C**** U component
3146.         DO 530 J=2,JM-1
3147.         DO 530 I=1,IM
3148.         FLUX = DT16*(FXY(I,J-1)+FXY(I,J))*(UA(I,J,L)+UA(I,J,L+1))
3149.         DUM(I,J,L)   = DUM(I,J,L)   - FLUX
3150.     530 DUM(I,J,L+1) = DUM(I,J,L+1) + FLUX
3151.         FLUX = DT16*FXYS*(UA(1,1,L)+UA(1,1,L+1))
3152.         DUM(1,1,L)   = DUM(1,1,L)   - FLUX
3153.         DUM(1,1,L+1) = DUM(1,1,L+1) + FLUX
3154.         FLUX = DT16*FXYN*(UA(1,JM,L)+UA(1,JM,L+1))
3155.         DUM(1,JM,L)  = DUM(1,JM,L)  - FLUX
3156.         DUM(1,JM,L+1)= DUM(1,JM,L+1)+ FLUX
3157.   C**** V component
3158.         IM1=IM
3159.         DO 550 J=1,JM-1
3160.         DO 550 I=1,IM
3161.         FLUX = DT16*(FXY(IM1,J)+FXY(I,J))*(VA(I,J,L)+VA(I,J,L+1))
3162.         DVM(I,J,L)   = DVM(I,J,L)   - FLUX
3163.         DVM(I,J,L+1) = DVM(I,J,L+1) + FLUX
3164.     550 IM1=I
3165.     560 continue
3166.   C****
3167.   C**** Add changes to momentum
3168.   C****
3169.         DO 640 L=1,LMA
3170.   C**** U component
3171.         DO 610 I=IM+1,IM*(JM-1)
3172.     610 UM2(I,1,L) = UM0(I,1,L) + DUM(I,1,L)
3173.         DUMS = 0.
3174.         DUMN = 0.
3175.         DO 620 I=1,IM
3176.         DUMS = DUMS + DUM(I, 1,L)
3177.     620 DUMN = DUMN + DUM(I,JM,L)
3178.         UM2(IM,1,L) = UM0(IM,1,L) + DUMS
3179.         UM2(1,JM,L) = UM0(1,JM,L) + DUMN
3180.   C**** V component
3181.         DO 630 I=1,IM*(JM-1)
3182.     630 VM2(I,1,L) = VM0(I,1,L) + DVM(I,1,L)
3183.     640 continue
3184.         RETURN
3185.         END
3500.   
3501.         SUBROUTINE APGF (UM2,VM2,H0M,HZM,MM,DT1)
3502.   C****
3503.   C**** APGF adds the Pressure Gradient Force to atmospheric momentum
3504.   C**** Input:  H0M (J), HZM (J), MM (kg), DT1 (s), MA (kg/m**2)
3505.   C**** Output: UM2 (kg*m/s) = UM2 - DT*(M*D(ZG+H*PK)-PK*D(H)))*DYP
3506.   C****         VM2 (kg*m/s) = VM2 - DT*(M*D(ZG+H*PK)-PK*D(H)))*DXV
3507.   C****
3508.         IMPLICIT REAL*8 (A-H,M,O-Z)
3509.         PARAMETER (IM=72,JM=46,LMA=9)
3510.         REAL*8 UM2(IM,JM,LMA),VM2(IM,JM,LMA),
3511.        *       H0M(IM,JM,LMA),HZM(IM,JM,LMA),MM(IM,JM,LMA)
3512.         COMMON /FIXDCB/ FTYPE(IM,JM,6),ZATMO(IM,JM)
3513.         COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM),DYV(JM)
3515.         COMMON /ATMOCB/ MA(IM,JM,LMA)
3515.5        COMMON /SURFCB/ BLDATA(IM,JM,8)
3516.         COMMON /WORK03/ DUM(IM,JM,LMA),DVM(IM,JM,LMA)
3517.         COMMON /WORK04/ PK(IM,JM,LMA),E(IM,JM,LMA),H(IM,JM,LMA),DZG(LMA)
3518.         PARAMETER (GRAV=9.81d0, RDRY=287., SHCD=1003.5, MSTRAT=101.9368d0,
3519.        *           RKAP=RDRY/SHCD, RKAPP1=RKAP+1., RKAPP2=RKAP+2.)
3520.   C****
3520.1  C**** FIXDCB  ZATMO   Atmospheric topography (m)
3520.2  C****
3520.3  C**** ATMOCB  MA      Atmospheric layer mass (kg/m**2)
3520.4  C****
3520.5  C**** SURFCB  BLDATA(5)  Atmospheric surface pressure (Pa-101325)
3520.6  C****         BLDATA(7)  Atmospheric surface geopotential (m^2/s^2)
3520.7  C****
3521.         DT2 = DT1/2.
3522.   C****
3523.   C**** Calculate the mass weighted geopotential ZG (m^2/s^2) plus
3524.   C**** sensible heat H*P^RKAP and the mass weighted P^RKAP
3525.   C****
3529.         DO 130 I=IM,IM*(JM-1)+1
3529.5  C**** Integrate pressures from the top down
3530.         MUP  = MSTRAT
3531.         PKUP = (MUP*GRAV)**RKAP
3534.         DO 110 L=LMA,1,-1
3535.         H(I,1,L) =  H0M(I,1,L)/MM(I,1,L)
3536.         HZ2    = 2.*HZM(I,1,L)/MM(I,1,L)
3537.         PKMUP  = PKUP*MUP
3538.         PKMMUP = PKMUP*MUP
3540.         M0     = MUP + .5*MA(I,1,L)
3541.         Y      = H(I,1,L) + HZ2*M0/MA(I,1,L)
3542.         MDN    = MUP + MA(I,1,L)
3543.         PKDN   = (MDN*GRAV)**RKAP
3544.         PKMDN  = PKDN*MDN
3545.         PKMMDN = PKMDN*MDN
3546.   C**** Mass weighted H*P^RKAP + ZG - ZGDN in the layer
3547.         E(I,1,L) = Y*PKDN - HZ2*(RKAP*PKMDN+
3548.        +  (PKMMDN-PKMMUP)/(MA(I,1,L)*RKAPP2))/(MA(I,1,L)*RKAPP1)
3549.   C**** Mass weighted P**KAPPA in the layer
3550.         PK(I,1,L) = (PKMDN-PKMUP)/(MA(I,1,L)*RKAPP1)
3551.   C**** Calaculate depth of ZG in the layer
3552.         DZG(L) = Y*(PKDN-PKUP) - HZ2*RKAP*PK(I,1,L)
3553.         MUP  = MDN
3554.     110 PKUP = PKDN
3554.1        BLDATA(I,1,5) = MDN*GRAV - 101325.
3554.2  C**** Integrate ZGDN from the bottom up, calculate H*P^RKAP + ZG
3554.3        ZGDN = ZATMO(I,1)*GRAV
3554.4  C     IF(FOCEAN(I,1).gt.0.)  ZGDN = BLDATA(I,1,7)
3554.5        E(I,1,1) = E(I,1,1) + ZGDN
3554.6        DO 120 L=2,LMA
3554.7        ZGDN = ZGDN + DZG(L-1)
3554.8    120 E(I,1,L) = E(I,1,L) + ZGDN
3555.     130 continue
3556.   C**** Replicate values at the poles
3557.         DO 140 L=1,LMA
3558.         DO 140 I=1,IM
3559.         E(I, 1,L) = E(IM,1,L)
3560.         E(I,JM,L) = E(1,JM,L)
3561.         H(I, 1,L) = H(IM,1,L)
3562.         H(I,JM,L) = H(1,JM,L)
3563.         PK(I,1,L) = PK(IM,1,L)
3564.     140 PK(I,JM,L)= PK(1,JM,L)
3565.   C****
3566.   C**** Calculate gradient of Pressure Force
3567.   C****
3568.         DO 440 L=1,LMA
3569.   C**** Smoothed East-West derivative affects the U-component
3570.         I=IM
3571.         DO 410 J=2,JM-1
3572.         DO 410 IP1=1,IM
3573.         DUM(I,J,L) =                -(E(IP1,J,L)-E(I,J,L) -
3574.        -  .5*(PK(IP1,J,L)+PK(I,J,L))*(H(IP1,J,L)-H(I,J,L)))*DT2*DYP(J)
3575.     410 I=IP1
3576.         CALL APFIL (DUM(1,1,L))
3577.   C     I=IM
3578.         DO 420 J=2,JM-1
3579.         DO 420 IP1=1,IM
3580.         DUM(I,J,L) = DUM(I,J,L)*(MA(IP1,J,L)+MA(I,J,L))
3581.     420 I=IP1
3582.   C**** North-South derivative affects the V-component
3583.         DO 430 J=1,JM-1
3584.         DO 430 I=1,IM
3585.     430 DVM(I,J,L) = -(MA(I,J+1,L)+MA(I,J,L))*(E(I,J+1,L)-E(I,J,L) -
3586.        -  .5*(PK(I,J+1,L)+PK(I,J,L))*(H(I,J+1,L)-H(I,J,L)))*DT2*DXV(J)
3587.     440 continue
3588.   C****
3589.   C**** Add pressure gradient force to momentum
3590.   C****
3591.         DO 940 L=1,LMA
3592.         DO 910 I=IM+1,IM*(JM-1)
3593.     910 UM2(I,1,L) = UM2(I,1,L) + DUM(I,1,L)
3594.         DO 930 I=1,IM*(JM-1)
3595.     930 VM2(I,1,L) = VM2(I,1,L) + DVM(I,1,L)
3596.     940 continue
3597.         RETURN
3598.         END
4000.   
4001.         SUBROUTINE AADVT (RM,RX,RY,RZ,MB,DT,QLIMIT,AIJL)
4002.   C****
4003.   C**** AADVT advects tracers using the linear upstream scheme.
4004.   C****
4005.   C**** Input:  MB (kg) = mass before advection
4006.   C****          DT (s) = time step
4007.   C****       MU (kg/s) = west to east mass flux
4008.   C****       MV (kg/s) = south to north mass flux
4009.   C****       MW (kg/s) = upward vertical mass flux
4010.   C****          QLIMIT = whether slope limitations should be used
4011.   C**** Output: RM (kg) = tracer mass
4012.   C****   RX,RY,RZ (kg) = first moments of tracer mass
4013.   C****       AIJL (kg) = diagnostic accumulation of tracer mass flux
4014.   C****
4015.         IMPLICIT REAL*8 (A-H,M,O-Z)
4016.         PARAMETER (IM=72,JM=46,LMA=9)
4017.         REAL*8 RM(IM,JM,LMA),RX(IM,JM,LMA),RY(IM,JM,LMA),RZ(IM,JM,LMA),
4018.        *       MB(IM,JM,LMA),AIJL(IM,JM,LMA,3)
4019.         COMMON /WORK03/ MA(IM,JM,LMA)
4020.   C****
4021.   C**** Load mass after advection from mass before advection
4022.   C****
4023.         DO 110 L=1,LMA
4024.         DO 110 I=IM,IM*(JM-1)+1
4025.     110 MA(I,1,L) = MB(I,1,L)
4026.   C****
4027.   C**** Advect the tracer using the Linear Upstream Scheme
4028.   C****
4029.     200 CALL AADVTX (RM,RX,RY,RZ,MA,.5*DT,QLIMIT,AIJL)
4030.         CALL AADVTY (RM,RX,RY,RZ,MA,   DT,QLIMIT,AIJL(1,1,1,2))
4031.         CALL AADVTZ (RM,RX,RY,RZ,MA,   DT,QLIMIT,AIJL(1,1,1,3))
4032.         CALL AADVTX (RM,RX,RY,RZ,MA,.5*DT,QLIMIT,AIJL)
4033.   C**** Fill in values at the poles
4034.         DO 210 L=1,LMA
4035.         DO 210 I=1,IM
4036.         RM(I, 1,L) = RM(IM,1,L)
4037.         RZ(I, 1,L) = RZ(IM,1,L)
4038.         RM(I,JM,L) = RM(1,JM,L)
4039.     210 RZ(I,JM,L) = RZ(1,JM,L)
4040.         RETURN
4041.         END
4200.   
4201.         SUBROUTINE AADVTX (RM,RX,RY,RZ,M,DT,QLIMIT,AIJL)
4202.   C****
4203.   C**** AADVTX advects tracers in the west to east direction using the
4204.   C**** linear upstream scheme.  If QLIMIT is true, the gradients are
4205.   C**** limited to prevent the mean tracer from becoming negative.
4206.   C****
4207.   C**** Input: DT (s) = time step
4208.   C****     MU (kg/s) = west to east mass flux
4209.   C****        QLIMIT = whether slope limitations should be used
4210.   C**** Input and Output: RM (kg) = tracer mass
4211.   C****             RX,RY,RZ (kg) = first moments of tracer mass
4212.   C****                    M (kg) = air mass
4213.   C****
4213.5        IMPLICIT REAL*8 (A-H,M,O-Z)
4214.         PARAMETER (IM=72,JM=46,LMA=9)
4215.         REAL*8 RM(IM,JM,LMA),RX(IM,JM,LMA),RY(IM,JM,LMA),RZ(IM,JM,LMA),
4216.        *        M(IM,JM,LMA),AIJL(IM,JM,LMA)
4217.         LOGICAL*4 QLIMIT
4218.         COMMON /FLUXCB/ MU(IM,JM,LMA)
4219.         COMMON /WORK04/ AM(IM),A(IM),FM(IM),FX(IM),FY(IM),FZ(IM)
4220.   C**** Loop over layers and latitudes
4221.         DO 320 L=1,LMA
4222.         DO 320 J=2,JM-1
4223.   C****
4224.   C**** Calculate FM (kg), FX (kg**2), FY (kg) and FZ (kg)
4225.   C****
4226.         I=IM
4227.         DO 120 IP1=1,IM
4228.         AM(I) = DT*MU(I,J,L)
4229.         IF(AM(I).LT.0.)  GO TO 110
4230.   C**** Air mass flux is positive
4231.         A(I)  = AM(I)/M(I,J,L)
4232.         IF(A(I).GT.1.)  WRITE (6,*) 'A>1:',I,J,L,A(I),M(I,J,L)
4233.         FM(I) = A(I)*(RM(I,J,L)+(1.-A(I))*RX(I,J,L))
4234.         FX(I) = AM(I)*(A(I)*A(I)*RX(I,J,L)-3.*FM(I))
4235.         FY(I) = A(I)*RY(I,J,L)
4236.         FZ(I) = A(I)*RZ(I,J,L)
4237.         GO TO 120
4238.   C**** Air mass flux is negative
4239.     110 A(I)  = AM(I)/M(IP1,J,L)
4240.         IF(A(I).LT.-1.)  WRITE (6,*) 'A<-1:',I,J,L,A(I),M(IP1,J,L)
4241.         FM(I) = A(I)*(RM(IP1,J,L)-(1.+A(I))*RX(IP1,J,L))
4242.         FX(I) = AM(I)*(A(I)*A(I)*RX(IP1,J,L)-3.*FM(I))
4243.         FY(I) = A(I)*RY(IP1,J,L)
4244.         FZ(I) = A(I)*RZ(IP1,J,L)
4245.     120 I=IP1
4246.   C****
4247.   C**** Modify the tracer moments so that the tracer mass in each
4248.   C**** division is non-negative
4249.   C****
4250.         IF(.NOT.QLIMIT)  GO TO 300
4251.         IM1=IM
4252.         DO 290 I=1,IM
4253.         IF(A(IM1).GE.0.)  GO TO 240
4254.   C**** Air is leaving through the left edge: 2 or 3 divisions
4255.         IF(FM(IM1).le.0.)  GO TO 210
4256.   C**** Left most division is negative, RML = -FM(I-1) < 0: Case 2 or 4
4257.         RX(I,J,L) = RM(I,J,L)/(1.+A(IM1))
4258.         FM(IM1) = 0.
4259.         FX(IM1) = AM(IM1)*A(IM1)*A(IM1)*RX(I,J,L)
4260.         IF(A(I).le.0.)  GO TO 290
4261.         FM(I) = A(I)*(RM(I,J,L)+(1.-A(I))*RX(I,J,L))
4262.         FX(I) = AM(I)*(A(I)*A(I)*RX(I,J,L)-3.*FM(I))
4263.         GO TO 290
4264.   C**** Left most division is non-negative, RML = -FM(I-() > 0:
4265.   C**** Case 1, 3 or 5
4266.     210 IF(A(I).le.0.)  GO TO 230
4267.   C**** Air is leaving through the right edge: 3 divisions
4268.         IF(FM(I).GE.0.)  GO TO 290
4269.   C**** Right most division is negative, RMR = FM(I) < 0: Case 3 or 5
4270.     220 RX(I,J,L) = -RM(I,J,L)/(1.-A(I))
4271.         FM(I) = 0.
4272.         FX(I) = AM(I)*A(I)*A(I)*RX(I,J,L)
4273.         FM(IM1) = A(IM1)*(RM(I,J,L)-(1.+A(IM1))*RX(I,J,L))
4274.         FX(IM1) = AM(IM1)*(A(IM1)*A(IM1)*RX(I,J,L)-3.*FM(IM1))
4275.         GO TO 290
4276.   C**** No air is leaving through the right edge: 2 divisions
4277.     230 IF(RM(I,J,L)+FM(IM1).GE.0.)  GO TO 290
4278.   C**** Right most division is negative, RMR = RM(I,J)+FM(I-1) < 0: Case 3
4279.         RX(I,J,L) = RM(I,J,L)/A(IM1)
4280.         FM(IM1) = -RM(I,J,L)
4281.         FX(IM1) = AM(IM1)*(A(IM1)+3.)*RM(I,J,L)
4282.         GO TO 290
4283.   C**** No air is leaving through the left edge: 1 or 2 divisions
4284.     240 IF(A(I).le.0.)  GO TO 290
4285.   C**** Air is leaving through the right edge: 2 divisions
4286.         IF(FM(I).GE.0.)  GO TO 250
4287.   C**** Right most division is negative, RMR = FM(I) < 0: Case 3
4288.         RX(I,J,L) = -RM(I,J,L)/(1.-A(I))
4289.         FM(I) = 0.
4290.         FX(I) = AM(I)*A(I)*A(I)*RX(I,J,L)
4291.         GO TO 290
4292.   C**** Right most division is non-negative, RMR = FM(I) > 0: Case 1 or 2
4293.     250 IF(RM(I,J,L)-FM(I).GE.0.)  GO TO 290
4294.   C**** Left most division is negative, RML = RM(I,J)-FM(I) < 0: Case 2
4295.         RX(I,J,L) = RM(I,J,L)/A(I)
4296.         FM(I) = RM(I,J,L)
4297.         FX(I) = AM(I)*(A(I)-3.)*RM(I,J,L)
4298.   C****
4299.     290 IM1=I
4300.   C****
4301.   C**** Calculate new tracer mass and first moments of tracer mass
4302.   C****
4303.     300 IM1=IM
4304.         DO 310 I=1,IM
4305.         RM(I,J,L) = RM(I,J,L) + (FM(IM1)-FM(I))
4306.         RX(I,J,L) = (RX(I,J,L)*M(I,J,L) + (FX(IM1)-FX(I))
4307.        *  + 3.*((AM(IM1)+AM(I))*RM(I,J,L)-M(I,J,L)*(FM(IM1)+FM(I))))
4308.        *  / (M(I,J,L)+(AM(IM1)-AM(I)))
4309.         RY(I,J,L) = RY(I,J,L) + (FY(IM1)-FY(I))
4310.         RZ(I,J,L) = RZ(I,J,L) + (FZ(IM1)-FZ(I))
4311.          M(I,J,L) =  M(I,J,L) + (AM(IM1)-AM(I))
4312.            IF(M(I,J,L).LT.0.)  GO TO 800
4313.            IF(QLIMIT .AND. RM(I,J,L).LT.0.)  GO TO 810
4314.            AIJL(I,J,L) = AIJL(I,J,L) + FM(I)
4315.     310 IM1=I
4316.     320 continue
4317.         RETURN
4318.     800 WRITE (6,*) 'M<0 in AADVTX:',I,J,L,M(I,J,L)
4319.     810 WRITE (6,*) 'RM in AADVTX:',I,J,L,RM(I,J,L)
4320.         WRITE (6,*) 'A=',(I,A(I),I=1,IM)
4321.         STOP
4322.         END
4400.   
4401.         SUBROUTINE AADVTY (RM,RX,RY,RZ,M,DT,QLIMIT,AIJL)
4402.   C****
4403.   C**** AADVTY advects tracers in the south to north direction using the
4404.   C**** linear upstream scheme.  If QLIMIT is true, the gradients are
4405.   C**** limited to prevent the mean tracer from becoming negative.
4406.   C****
4407.   C**** Input: DT (s) = time step
4408.   C****     MV (kg/s) = south to north mass flux
4409.   C****        QLIMIT = whether slope limitations should be used
4410.   C**** Input and Output: RM (kg) = tracer mass
4411.   C****             RX,RY,RZ (kg) = first moments of tracer mass
4412.   C****                    M (kg) = air mass
4413.   C****
4413.5        IMPLICIT REAL*8 (A-H,M,O-Z)
4414.         PARAMETER (IM=72,JM=46,LMA=9)
4415.         REAL*8 RM(IM,JM,LMA),RX(IM,JM,LMA),RY(IM,JM,LMA),RZ(IM,JM,LMA),
4416.        *        M(IM,JM,LMA),AIJL(IM,JM,LMA)
4417.         LOGICAL*4 QLIMIT
4418.         COMMON /FLUXCB/ MU(IM,JM,LMA),MV(IM,JM,LMA)
4419.         COMMON /WORK04/ BM(JM),B(JM),FM(JM),FX(JM),FY(JM),FZ(JM)
4420.   C**** Loop over layers and longitudes
4421.         DO 340 L=1,LMA
4422.         SBMS = 0.
4423.         SFMS = 0.
4424.         SFZS = 0.
4425.         SBMN = 0.
4426.         SFMN = 0.
4427.         SFZN = 0.
4428.         DO 330 I=1,IM
4429.   C****
4430.   C**** Calculate FM (kg), FX (kg), FY (kg**2) and FZ (kg) near South Pole
4431.   C****
4432.         BM(1) = DT*MV(I,1,L)
4433.         IF(BM(1).LT.0.)  GO TO 110
4434.   C**** Air mass flux is positive
4435.         B(1)  = BM(1)/M(IM,1,L)
4436.         IF(B(1).GT.1.)  WRITE (6,*) 'B>1:',I,1,L,B(1),M(IM,1,L)
4437.         FM(1) = B(1)*RM(IM,1,L)
4438.         FX(1) = 0.
4439.         FY(1) = -3.*BM(1)*FM(1)
4440.         FZ(1) = B(1)*RZ(IM,1,L)
4441.         GO TO 120
4442.   C**** Air mass flux is negative
4443.     110 B(1)  = BM(1)/M(I,2,L)
4444.         IF(B(1).LT.-1.)  WRITE (6,*) 'B<-1:',I,1,L,B(1),M(1,2,L)
4445.         FM(1) = B(1)*(RM(I,2,L)-(1.+B(1))*RY(I,2,L))
4446.         FX(1) = B(1)*RX(I,2,L)
4447.         FY(1) = BM(1)*(B(1)*B(1)*RY(I,2,L)-3.*FM(1))
4448.         FZ(1) = B(1)*RZ(I,2,L)
4449.   C****
4450.   C**** Calculate FM (kg), FX (kg), FY (kg**2) and FZ (kg) in the interior
4451.   C****
4452.     120 DO 140 J=2,JM-2
4453.         BM(J) = DT*MV(I,J,L)
4454.         IF(BM(J).LT.0.)  GO TO 130
4455.   C**** Air mass flux is positive
4456.         B(J)  = BM(J)/M(I,J,L)
4457.         IF(B(J).GT.1.)  WRITE (6,*) 'B>1:',I,J,L,B(J),M(I,J,L)
4458.         FM(J) = B(J)*(RM(I,J,L)+(1.-B(J))*RY(I,J,L))
4459.         FX(J) = B(J)*RX(I,J,L)
4460.         FY(J) = BM(J)*(B(J)*B(J)*RY(I,J,L)-3.*FM(J))
4461.         FZ(J) = B(J)*RZ(I,J,L)
4462.         GO TO 140
4463.   C**** Air mass flux is negative
4464.     130 B(J)  = BM(J)/M(I,J+1,L)
4465.         IF(B(J).LT.-1.)  WRITE (6,*) 'B<-1:',I,J,L,B(J),M(I,J+1,L)
4466.         FM(J) = B(J)*(RM(I,J+1,L)-(1.+B(J))*RY(I,J+1,L))
4467.         FX(J) = B(J)*RX(I,J+1,L)
4468.         FY(J) = BM(J)*(B(J)*B(J)*RY(I,J+1,L)-3.*FM(J))
4469.         FZ(J) = B(J)*RZ(I,J+1,L)
4470.     140 continue
4471.   C****
4472.   C**** Calculate FM (kg), FX (kg), FY (kg**2) and FZ (kg) near North Pole
4473.   C****
4474.         BM(JM-1) = DT*MV(I,JM-1,L)
4475.         IF(BM(JM-1).LT.0.)  GO TO 150
4476.   C**** Air mass flux is positive
4477.         B(JM-1)  = BM(JM-1)/M(I,JM-1,L)
4478.         IF(B(JM-1).GT.1.)  WRITE(6,*) 'B>1:',I,JM-1,L,B(JM-1),M(I,JM-1,L)
4479.         FM(JM-1) = B(JM-1)*(RM(I,JM-1,L)+(1.-B(JM-1))*RY(I,JM-1,L))
4480.         FX(JM-1) = B(JM-1)*RX(I,JM-1,L)
4481.         FY(JM-1) = BM(JM-1)*(B(JM-1)*B(JM-1)*RY(I,JM-1,L)-3.*FM(JM-1))
4482.         FZ(JM-1) = B(JM-1)*RZ(I,JM-1,L)
4483.         GO TO 200
4484.   C**** Air mass flux is negative
4485.     150 B(JM-1)  = BM(JM-1)/M(1,JM,L)
4486.         IF(B(JM-1).LT.-1.)  WRITE (6,*) 'B<-1:',I,JM-1,L,B(JM-1),M(1,JM,L)
4487.         FM(JM-1) = B(JM-1)*RM(1,JM,L)
4488.         FX(JM-1) = 0.
4489.         FY(JM-1) = -3.*BM(JM-1)*FM(JM-1)
4490.         FZ(JM-1) = B(JM-1)*RZ(1,JM,L)
4491.   C****
4492.   C**** Modify the tracer moments so that the tracer mass in each
4493.   C**** division is non-negative
4494.   C****
4495.     200 IF(.NOT.QLIMIT)  GO TO 300
4496.         DO 290 J=2,JM-1
4497.         IF(B(J-1).GE.0.)  GO TO 240
4498.   C**** Air is leaving through the south edge: 2 or 3 divisions
4499.         IF(FM(J-1).le.0.)  GO TO 210
4500.   C**** South most division is negative, RMS = -FM(J-1) < 0: Case 2 or 4
4501.         RY(I,J,L) = RM(I,J,L)/(1.+B(J-1))
4502.         FM(J-1) = 0.
4503.         FY(J-1) = BM(J-1)*B(J-1)*B(J-1)*RY(I,J,L)
4504.         IF(B(J).le.0.)  GO TO 290
4505.         FM(J) = B(J)*(RM(I,J,L)+(1.-B(J))*RY(I,J,L))
4506.         FY(J) = BM(J)*(B(J)*B(J)*RY(I,J,L)-3.*FM(J))
4507.         GO TO 290
4508.   C**** South most division is non-negative, RMS = -FM(J-1) > 0:
4509.   C**** Case 1, 3 or 5
4510.     210 IF(B(J).le.0.)  GO TO 230
4511.   C**** Air is leaving through the north edge: 3 divisions
4512.         IF(FM(J).GE.0.)  GO TO 290
4513.   C**** North most division is negative, RMN = FM(J) < 0: Case 3 or 5
4514.     220 RY(I,J,L) = -RM(I,J,L)/(1.-B(J))
4515.         FM(J) = 0.
4516.         FY(J) = BM(J)*B(J)*B(J)*RY(I,J,L)
4517.         FM(J-1) = B(J-1)*(RM(I,J,L)-(1.+B(J-1))*RY(I,J,L))
4518.         FY(J-1) = BM(J-1)*(B(J-1)*B(J-1)*RY(I,J,L)-3.*FM(J-1))
4519.         GO TO 290
4520.   C**** No air is leaving through the north edge: 2 divisions
4521.     230 IF(RM(I,J,L)+FM(J-1).GE.0.)  GO TO 290
4522.   C**** North most division is negative, RMN = RM(I,J)+FM(J-1) < 0: Case 3
4523.         RY(I,J,L) = RM(I,J,L)/B(J-1)
4524.         FM(J-1) = -RM(I,J,L)
4525.         FY(J-1) = BM(J-1)*(B(J-1)+3.)*RM(I,J,L)
4526.         GO TO 290
4527.   C**** No air is leaving through the south edge: 1 or 2 divisions
4528.     240 IF(B(J).le.0.)  GO TO 290
4529.   C**** Air is leaving through the north edge: 2 divisions
4530.         IF(FM(J).GE.0.)  GO TO 250
4531.   C**** North most division is negative, RMN = FM(J) < 0: Case 3
4532.         RY(I,J,L) = -RM(I,J,L)/(1.-B(J))
4533.         FM(J) = 0.
4534.         FY(J) = BM(J)*B(J)*B(J)*RY(I,J,L)
4535.         GO TO 290
4536.   C**** North most division is non-negative, RMN = FM(J) > 0: Case 1 or 2
4537.     250 IF(RM(I,J,L)-FM(J).GE.0.)  GO TO 290
4538.   C**** South most division is negative, RMS = RM(I,J)-FM(J) < 0: Case 2
4539.         RY(I,J,L) = RM(I,J,L)/B(J)
4540.         FM(J) = RM(I,J,L)
4541.         FY(J) = BM(J)*(B(J)-3.)*RM(I,J,L)
4542.   C****
4543.     290 continue
4544.     300 SBMS = SBMS + BM(1)
4545.         SFMS = SFMS + FM(1)
4546.         SFZS = SFZS + FZ(1)
4547.         SBMN = SBMN + BM(JM-1)
4548.         SFMN = SFMN + FM(JM-1)
4549.         SFZN = SFZN + FZ(JM-1)
4550.   C****
4551.   C**** Calculate new tracer mass and first moments of tracer mass
4552.   C****
4553.         DO 310 J=2,JM-1
4554.         RM(I,J,L) = RM(I,J,L) + (FM(J-1)-FM(J))
4555.         RX(I,J,L) = RX(I,J,L) + (FX(J-1)-FX(J))
4556.         RY(I,J,L) = (RY(I,J,L)*M(I,J,L) + (FY(J-1)-FY(J))
4557.        *  + 3.*((BM(J-1)+BM(J))*RM(I,J,L)-M(I,J,L)*(FM(J-1)+FM(J))))
4558.        *  / (M(I,J,L)+(BM(J-1)-BM(J)))
4559.         RZ(I,J,L) = RZ(I,J,L) + (FZ(J-1)-FZ(J))
4560.          M(I,J,L) =  M(I,J,L) + (BM(J-1)-BM(J))
4561.            IF(M(I,J,L).le.0.)  GO TO 800
4562.            IF(QLIMIT .AND. RM(I,J,L).LT.0.)  GO TO 810
4563.     310 continue
4564.            DO 320 J=1,JM-1
4565.     320    AIJL(I,J,L) = AIJL(I,J,L) + FM(J)
4566.     330    continue
4567.         RM(IM,1,L) = RM(IM,1,L) - SFMS/IM
4568.         RZ(IM,1,L) = RZ(IM,1,L) - SFZS/IM
4569.          M(IM,1,L) =  M(IM,1,L) - SBMS/IM
4570.         RM(1,JM,L) = RM(1,JM,L) + SFMN/IM
4571.         RZ(1,JM,L) = RZ(1,JM,L) + SFZN/IM
4572.          M(1,JM,L) =  M(1,JM,L) + SBMN/IM
4573.            IF(M(IM,1,L).le.0.)  GO TO 800
4574.            IF(M(1,JM,L).le.0.)  GO TO 800
4575.            IF(QLIMIT .AND. RM(IM,1,L).LT.0.)  GO TO 810
4576.            IF(QLIMIT .AND. RM(1,JM,L).LT.0.)  GO TO 810
4577.     340 continue
4578.         RETURN
4579.     800 WRITE (6,*) 'M<0 in AADVTY:',I,J,L,M(I,J,L)
4580.     810 WRITE (6,*) 'RM in AADVTY:',I,J,L,RM(I,J,L)
4581.         WRITE (6,*) 'B=',(J,B(J),J=1,JM-1)
4582.         STOP
4583.         END
4600.   
4601.         SUBROUTINE AADVTZ (RM,RX,RY,RZ,M,DT,QLIMIT,AIJL)
4602.   C****
4603.   C**** AADVTZ advects tracers in the vertical direction using the
4604.   C**** linear upstream scheme.  If QLIMIT is true, the gradients are
4605.   C**** limited to prevent the mean tracer from becoming negative.
4606.   C****
4607.   C**** Input: DT (s) = time step
4608.   C****     MW (kg/s) = vertical mass flux
4609.   C****        QLIMIT = whether slope limitations should be used
4610.   C**** Input and Output: RM (kg) = tracer mass
4611.   C****             RX,RY,RZ (kg) = first moments of tracer mass
4612.   C****                    M (kg) = air mass
4613.   C****
4613.5        IMPLICIT REAL*8 (A-H,M,O-Z)
4614.         PARAMETER (IM=72,JM=46,LMA=9)
4615.         REAL*8 RM(IM,JM,LMA),RX(IM,JM,LMA),RY(IM,JM,LMA),RZ(IM,JM,LMA),
4616.        *        M(IM,JM,LMA),AIJL(IM,JM,LMA)
4617.         LOGICAL*4 QLIMIT
4618.         COMMON /FLUXCB/ MU(IM,JM,LMA),MV(IM,JM,LMA),MW(IM,JM,LMA-1)
4619.         COMMON /WORK04/ CM(LMA),C(0:LMA),
4620.        *  FM(0:LMA),FX(LMA),FY(LMA),FZ(0:LMA)
4621.   C**** Loop over latitudes and longitudes
4622.         DO 330 I=IM,IM*(JM-1)+1
4623.   C****
4624.   C**** Calculate FM (kg), FX (kg), FY (kg) and FZ (kg**2)
4625.   C****
4626.         DO 120 L=1,LMA-1
4627.         CM(L) = DT*MW(I,1,L)
4628.         IF(CM(L).LT.0.) GO TO 110
4629.   C**** Air mass flux is positive
4630.         C(L)  = CM(L)/M(I,1,L)
4631.         IF(C(L).GT.1.)  WRITE (6,*) 'C>1:',I,L,C(L),M(I,1,L)
4632.         FM(L) = C(L)*(RM(I,1,L)+(1.-C(L))*RZ(I,1,L))
4633.         FX(L) = C(L)*RX(I,1,L)
4634.         FY(L) = C(L)*RY(I,1,L)
4635.         FZ(L) = CM(L)*(C(L)*C(L)*RZ(I,1,L)-3.*FM(L))
4636.         GO TO 120
4637.   C**** Air mass flux is negative
4638.     110 C(L)  = CM(L)/M(I,1,L+1)
4639.         IF(C(L).LT.-1.)  WRITE (6,*) 'C<-1:',I,L,C(L),M(I,1,L+1)
4640.         FM(L) = C(L)*(RM(I,1,L+1)-(1.+C(L))*RZ(I,1,L+1))
4641.         FX(L) = C(L)*RX(I,1,L+1)
4642.         FY(L) = C(L)*RY(I,1,L+1)
4643.         FZ(L) = CM(L)*(C(L)*C(L)*RZ(I,1,L+1)-3.*FM(L))
4644.     120 continue
4645.   C****
4646.   C**** Modify the tracer moments so that the tracer mass in each
4647.   C**** division is non-negative
4648.   C****
4649.         IF(.NOT.QLIMIT)  GO TO 300
4650.         C(0)  = 0.
4651.         C(LMA) = 0.
4652.         DO 290 L=1,LMA
4653.         IF(C(L-1).GE.0.)  GO TO 240
4654.   C**** Air is leaving through the bottom edge: 2 or 3 divisions
4655.         IF(FM(L-1).le.0.)  GO TO 210
4656.   C**** Bottom most division is negative, RMB = -FM(L-1) < 0: Case 2 or 4
4657.         RZ(I,1,L) = RM(I,1,L)/(1.+C(L-1))
4658.         FM(L-1) = 0.
4659.         FZ(L-1) = CM(L-1)*C(L-1)*C(L-1)*RZ(I,1,L)
4660.         IF(C(L).le.0.)  GO TO 290
4661.         FM(L) = C(L)*(RM(I,1,L)+(1.-C(L))*RZ(I,1,L))
4662.         FZ(L) = CM(L)*(C(L)*C(L)*RZ(I,1,L)-3.*FM(L))
4663.         GO TO 290
4664.   C**** Bottom most division is non-negative, RMB = -FM(L-1) > 0:
4665.   C**** Case 1, 3 or 5
4666.     210 IF(C(L).le.0.)  GO TO 230
4667.   C**** Air is leaving through the top edge: 3 divisions
4668.         IF(FM(L).GE.0.)  GO TO 290
4669.   C**** Top most division is negative, RMT = FM(L) < 0: Case 3 or 5
4670.     220 RZ(I,1,L) = -RM(I,1,L)/(1.-C(L))
4671.         FM(L) = 0.
4672.         FZ(L) = CM(L)*C(L)*C(L)*RZ(I,1,L)
4673.         FM(L-1) = C(L-1)*(RM(I,1,L)-(1.+C(L-1))*RZ(I,1,L))
4674.         FZ(L-1) = CM(L-1)*(C(L-1)*C(L-1)*RZ(I,1,L)-3.*FM(L-1))
4675.         GO TO 290
4676.   C**** No air is leaving through the top edge: 2 divisions
4677.     230 IF(RM(I,1,L)+FM(L-1).GE.0.)  GO TO 290
4678.   C**** Top most division is negative, RMT = RM(I,1,L)+FM(L-1) < 0: Case 3
4679.         RZ(I,1,L) = RM(I,1,L)/C(L-1)
4680.         FM(L-1) = -RM(I,1,L)
4681.         FZ(L-1) = CM(L-1)*(C(L-1)+3.)*RM(I,1,L)
4682.         GO TO 290
4683.   C**** No air is leaving through the bottom edge: 1 or 2 divisions
4684.     240 IF(C(L).le.0.)  GO TO 290
4685.   C**** Air is leaving through the top edge: 2 divisions
4686.         IF(FM(L).GE.0.)  GO TO 250
4687.   C**** Top most division is negative, RMT = FM(L) < 0: Case 3
4688.         RZ(I,1,L) = -RM(I,1,L)/(1.-C(L))
4689.         FM(L) = 0.
4690.         FZ(L) = CM(L)*C(L)*C(L)*RZ(I,1,L)
4691.         GO TO 290
4692.   C**** Top most division is non-negative, RMT = FM(L) > 0: Case 1 or 2
4693.     250 IF(RM(I,1,L)-FM(L).GE.0.)  GO TO 290
4694.   C**** Bottom most division is negative, RMB = RM(I,1,L)-FM(L) < 0: Cas 2
4695.         RZ(I,1,L) = RM(I,1,L)/C(L)
4696.         FM(L) = RM(I,1,L)
4697.         FZ(L) = CM(L)*(C(L)-3.)*RM(I,1,L)
4698.   C****
4699.     290 continue
4700.   C****
4701.   C**** Calculate new tracer mass and first moments of tracer mass
4702.   C****
4703.   C**** Calculation in the first layer
4704.     300 RM(I,1,1) = RM(I,1,1) - FM(1)
4705.         RX(I,1,1) = RX(I,1,1) - FX(1)
4706.         RY(I,1,1) = RY(I,1,1) - FY(1)
4707.         RZ(I,1,1) = (RZ(I,1,1)*M(I,1,1) - FZ(1)
4708.        *  + 3.*(CM(1)*RM(I,1,1)-M(I,1,1)*FM(1)))
4709.        *  / (M(I,1,1)-CM(1))
4710.          M(I,1,1) =  M(I,1,1) - CM(1)
4711.            IF(M(I,1,1).le.0.)  GO TO 800
4712.            IF(QLIMIT .AND. RM(I,1,1).LT.0.)  GO TO 810
4713.   C**** Calculation in the interior layers
4714.         DO 310 L=2,LMA-1
4715.         RM(I,1,L) = RM(I,1,L) + (FM(L-1)-FM(L))
4716.         RX(I,1,L) = RX(I,1,L) + (FX(L-1)-FX(L))
4717.         RY(I,1,L) = RY(I,1,L) + (FY(L-1)-FY(L))
4718.         RZ(I,1,L) = (RZ(I,1,L)*M(I,1,L) + (FZ(L-1)-FZ(L))
4719.        *  + 3.*((CM(L-1)+CM(L))*RM(I,1,L)-M(I,1,L)*(FM(L-1)+FM(L))))
4720.        *  / (M(I,1,L)+(CM(L-1)-CM(L)))
4721.          M(I,1,L) =  M(I,1,L) + (CM(L-1)-CM(L))
4722.            IF(M(I,1,L).le.0.)  GO TO 800
4723.            IF(QLIMIT .AND. RM(I,1,L).LT.0.)  GO TO 810
4724.     310 continue
4725.   C**** Calculation in the top layer
4726.         RM(I,1,LMA) = RM(I,1,LMA) + FM(LMA-1)
4727.         RX(I,1,LMA) = RX(I,1,LMA) + FX(LMA-1)
4728.         RY(I,1,LMA) = RY(I,1,LMA) + FY(LMA-1)
4729.         RZ(I,1,LMA) = (RZ(I,1,LMA)*M(I,1,LMA) + FZ(LMA-1)
4730.        *  + 3.*(CM(LMA-1)*RM(I,1,LMA)-M(I,1,LMA)*FM(LMA-1)))
4731.        *  / (M(I,1,LMA)+CM(LMA-1))
4732.          M(I,1,LMA) =  M(I,1,LMA) + CM(LMA-1)
4733.            IF(M(I,1,LMA).le.0.)  GO TO 800
4734.            IF(QLIMIT .AND. RM(I,1,LMA).LT.0.)  GO TO 810
4735.            DO 320 L=1,LMA-1
4736.     320    AIJL(I,1,L) = AIJL(I,1,L) + FM(L)
4737.     330 continue
4738.         RETURN
4739.     800 WRITE (6,*) 'M<0 in AADVTZ:',I,L,M(I,1,L)
4740.     810 WRITE (6,*) 'RM in AADVTZ:',I,L,RM(I,1,L)
4741.         WRITE (6,*) 'C=',(L,C(L),L=0,LMA)
4742.         STOP
4743.         END
5000.   
5001.         SUBROUTINE AABFIL
5002.   C****
5003.   C**** AABFIL applies an 8-th order Alternating Binomial FILter to
5004.   C**** Atmospheric fields UA and VA.  The frequency and strength depends
5005.   C**** on NFILXA and NFILYA.
5006.   C****
5007.         INCLUDE 'C070.COM'
5010.         LOGICAL*4 QFILX,QFILY
5011.         COMMON /WORK02/ MXY(IM,JM,LMA),MXYY(IM,JM,LMA),MXXY(IM,JM,LMA),
5011.1       *  X(IM),Y(JM)
5012.   C****
5013.   C**** Filtering in east-west direction
5014.   C****
5015.         QFILX = .FALSE.
5016.         IF(NFILXA.le.0)  GO TO 400
5017.         IF(MOD(IDT-IDT0,NFILXA).ne.0)  GO TO 400
5018.         QFILX = .TRUE.
5019.         X4TO8 = 1./(NFILXA*4.**8)
5019.5        DO 350 L=1,LMA
5020.   C**** Calculate mass weighting for UA and VA
5021.         DO 110 J=2,JM-1
5022.         DO 110 I=1,IM
5023.     110 MXYY(I,J,L) = MXY(I,J-1,L)+MXY(I,J,L)
5024.         IM1=IM
5025.         DO 120 J=1,JM-1
5026.         DO 120 I=1,IM
5027.         MXXY(I,J,L) = MXY(IM1,J,L)+MXY(I,J,L)
5028.     120 IM1=I
5029.   C**** Filter U component of momentum
5031.         DO 240 J=2,JM-1
5032.         DO 210 I=1,IM
5033.     210 X(I) = UA(I,J,L)*MXYY(I,J,L)
5034.         DO 230 N=1,8
5035.         X1   = X(1)
5036.         XIM1 = X(IM)
5037.         DO 220 I=1,IM-1
5038.         XI   = X(I)
5039.         X(I) = XIM1-XI-XI+X(I+1)
5040.     220 XIM1 = XI
5041.     230 X(IM)= XIM1-X(IM)-X(IM)+X1
5042.         DO 240 I=1,IM
5043.     240 UA(I,J,L) = UA(I,J,L) - X(I)*X4TO8/MXYY(I,J,L)
5044.   C**** Filter V component of momentum
5045.         DO 340 J=1,JM-1
5046.         DO 310 I=1,IM
5047.     310 X(I) = VA(I,J,L)*MXXY(I,J,L)
5048.         DO 330 N=1,8
5049.         X1   = X(1)
5050.         XIM1 = X(IM)
5051.         DO 320 I=1,IM-1
5052.         XI   = X(I)
5053.         X(I) = XIM1-XI-XI+X(I+1)
5054.     320 XIM1 = XI
5055.     330 X(IM)= XIM1-X(IM)-X(IM)+X1
5056.         DO 340 I=1,IM
5057.     340 VA(I,J,L) = VA(I,J,L) - X(I)*X4TO8/MXXY(I,J,L)
5058.     350 continue
5059.   C****
5060.   C**** Filtering in north-south direction
5061.   C****
5062.     400 QFILY = .FALSE.
5063.         IF(NFILYA.le.0)  GO TO 700
5064.         IF(MOD(IDT-IDT0,NFILYA).ne.0)  GO TO 700
5065.         QFILY = .TRUE.
5066.         Y4TO8 = 1./(NFILYA*4.**8)
5067.   C**** Calculate mass weighting for UA and VA
5067.5        DO 650 L=1,LMA
5068.         IF(QFILX)  GO TO 430
5069.         DO 410 J=2,JM-1
5070.         DO 410 I=1,IM
5071.     410 MXYY(I,J,L) = MXY(I,J-1,L)+MXY(I,J,L)
5072.         IM1=IM
5073.         DO 420 J=1,JM-1
5074.         DO 420 I=1,IM
5075.         MXXY(I,J,L) = MXY(IM1,J,L)+MXY(I,J,L)
5076.     420 IM1=I
5077.     430 DO 440 I=1,IM
5078.         MXYY(I, 1,L) = MXY(I, 1  ,L)
5079.     440 MXYY(I,JM,L) = MXY(I,JM-1,L)
5080.   C**** Filter U component of momentum
5082.         DO 540 I=1,IM
5083.         DO 510 J=1,JM
5084.     510 Y(J) = UA(I,J,L)*MXYY(I,J,L)
5085.         DO 530 N=1,8
5086.         YJM1 = Y(1)
5087.         Y(1) = Y(2)-Y(1)
5088.         DO 520 J=2,JM-1
5089.         YJ   = Y(J)
5090.         Y(J) = YJM1-YJ-YJ+Y(J+1)
5091.     520 YJM1 = YJ
5092.     530 Y(JM)= YJM1-Y(JM)
5093.         DO 540 J=1,JM
5094.     540 UA(I,J,L) = UA(I,J,L) - Y(J)*Y4TO8/MXYY(I,J,L)
5095.   C**** Make UA winds at poles to be uniform
5096.         DO 560 J=1,JM,JM-1
5097.         YJ = 0.
5098.         DO 550 I=1,IM
5099.     550 YJ = YJ + UA(I,J,L)
5100.         DO 560 I=1,IM
5101.     560 UA(I,J,L) = YJ/IM
5102.   C**** Filter V component of momentum
5103.         DO 640 I=1,IM
5104.         DO 610 J=1,JM-1
5105.     610 Y(J) = VA(I,J,L)*MXXY(I,J,L)
5106.         DO 630 N=1,8
5107.         YJM1 = Y(1)
5108.         Y(1) = Y(2)-Y(1)
5109.         DO 620 J=2,JM-2
5110.         YJ   = Y(J)
5111.         Y(J) = YJM1-YJ-YJ+Y(J+1)
5112.     620 YJM1 = YJ
5113.     630 Y(JM-1)= YJM1-Y(JM-1)
5114.         DO 640 J=1,JM-1
5115.     640 VA(I,J,L) = VA(I,J,L) - Y(J)*Y4TO8/MXXY(I,J,L)
5116.     650 continue
5117.   C****
5118.     700 RETURN
5119.         END
6000.   
6001.         SUBROUTINE DAILY
6002.   C****
6003.   C**** DAILY performs the following functions at the beginning
6004.   C**** of each new day:
6005.   C**** 1. Restore global air mass
6006.   C**** 2. Update daily calendar
6007.   C**** 3. Add daily minimum and maximum temperatures to AIJ arrays
6008.   C****
6009.         INCLUDE 'C070.COM'
6010.         PARAMETER (AREAG=2.*TWOPI*RADIUS*RADIUS)
6011.         INTEGER*4 JDOFM(13)
6012.         CHARACTER AMONTH(12)*4
6013.         DATA AMONTH/'Jan','Feb','Mar','Apr','May','Jun',
6014.        *            'Jul','Aug','Sep','Oct','Nov','Dec'/
6015.         DATA JDOFM /0,31,59,90,120,151,181,212,243,273,304,334,365/
6016.   C****
6017.   C**** Global air mass is kept constant at MDRYA (kg/m^2)
6018.   C****
6019.   C**** Calculate current global air mass
6020.     200 SMASS  = 0.
6021.         DO 220 L=1,LMA
6022.         DO 220 J=1,JM
6023.         SMASSJ = 0.
6024.         DO 210 I=1,IM
6025.     210 SMASSJ = SMASSJ + MA(I,J,L)
6026.     220 SMASS  = SMASS  + SMASSJ*DXYP(J)
6027.   C**** Correct air mass field
6028.         DELTAM = MDRYA - SMASS/AREAG - MSTRAT
6029.         DO 230 L=1,LMA
6030.         DO 230 I=1,IM*JM
6031.     230 MA(I,1,L) = MA(I,1,L) + DELTAM*DSIGA(L)
6032.         QPK  = .FALSE.
6033.         WRITE (6,923) DELTAM
6034.   C****
6035.   C**** Calculate the daily calendar
6036.   C****
6037.     300 JYEAR = IYEAR
6038.         JDAY  = IDAY+1 - (IYEAR-IYEAR0)*365
6039.         DO 310 KMONTH=1,12
6040.         IF(JDAY.le.JDOFM(KMONTH+1))  GO TO 320
6041.     310 continue
6042.     320 JDATE = JDAY-JDOFM(KMONTH)
6043.         JMON  = AMONTH(KMONTH)
6044.         IF(IDT.le.IDT0)  RETURN
6045.   C****
6046.   C**** Add minimum and maximum temperatures and daily variances to AIJ
6047.   C****
6049.         IDACC(12) = IDACC(12) + 1
6050.         AIJ(1,1,10) = AIJ(1,1,10) + TMNMX(1,1,1)
6051.         AIJ(1,1,11) = AIJ(1,1,11) + TMNMX(1,1,2)
6052.         AIJ(1,1,12) = AIJ(1,1,12) + TMNMX(1,1,3)
6053.         AIJ(1,1,13) = AIJ(1,1,13) + TMNMX(1,1,4)
6054.         AIJ(1,1,73) = AIJ(1,1,73) + TMNMX(1,1,5)*TMNMX(1,1,5)
6055.         AIJ(1,1,74) = AIJ(1,1,74) + TMNMX(1,1,6)*TMNMX(1,1,6)
6056.         AIJ(1,1,75) = AIJ(1,1,75) + TMNMX(IM,1,7)*TMNMX(IM,1,7)
6058.         DO 410 I=IM+1,IM*(JM-1)+1
6059.         AIJ(I,1,10) = AIJ(I,1,10) + TMNMX(I,1,1)
6060.         AIJ(I,1,11) = AIJ(I,1,11) + TMNMX(I,1,2)
6061.         AIJ(I,1,12) = AIJ(I,1,12) + TMNMX(I,1,3)
6062.         AIJ(I,1,13) = AIJ(I,1,13) + TMNMX(I,1,4)
6063.         AIJ(I,1,73) = AIJ(I,1,73) + TMNMX(I,1,5)*TMNMX(I,1,5)
6064.         AIJ(I,1,74) = AIJ(I,1,74) + TMNMX(I,1,6)*TMNMX(I,1,6)
6065.         AIJ(I,1,75) = AIJ(I,1,75) + TMNMX(I,1,7)*TMNMX(I,1,7)
6067.     410 continue
6068.         RETURN
6069.   C****
6070.     923 FORMAT ('0Air mass added in DAILY is:',F10.6/)
6071.         END
7000.   
7001.         SUBROUTINE CHECKT (SUBRN)
7002.   C****
7003.   C**** CHECKT checks whether the subsurface temperatures are reasonable.
7004.   C**** CHECKT is turned on in the run deck by setting QCHECK=T.
7005.   C****
7006.         INCLUDE 'C070.COM'
7007.         CHARACTER*6 SUBRN
7008.         IF(.not.QCHECK)  RETURN
7009.   C****
7010.   C**** Check ground and ice reservoirs for reasonable values
7011.   C****
7012.         DO 150 J=1,JM
7013.         IMAX=IM
7014.         IF(J.eq.1 .or. J.eq.JM)  IMAX=1
7015.         DO 150 I=1,IMAX
7016.         IF(FWATER(I,J).le.0. .or. RSI(I,J).le.0.)  GO TO 110
7017.   C**** Check temperatures over Sea Ice or Lake Ice
7018.         TSI1 = (HSI(I,J,1)/XSI1/MSI(I,J,1) + ELHM) / SHCI
7019.         TSI2 = (HSI(I,J,2)/XSI2/MSI(I,J,1) + ELHM) / SHCI
7020.         TSI3 = (HSI(I,J,3)/XSI3/MSI(I,J,2) + ELHM) / SHCI
7021.         TSI4 = (HSI(I,J,4)/XSI4/MSI(I,J,2) + ELHM) / SHCI
7022.         IF(TSI1.le.1.D-12 .and. TSI2.le.1.D-12 .and. TSI1.ge.-80. .and.
7023.        *   TSI3.le.1.D-12 .and. TSI4.le.1.D-12)  GO TO 110
7024.         WRITE (6,*) 'After ',SUBRN,': I,J,MSI,TSI=',I,J,
7025.        *      MSI(I,J,1),MSI(I,J,2),TSI1,TSI2,TSI3,TSI4
7026.         STOP 'CHECKT'
7027.     110 IF(FGICE(I,J).le.0.)  GO TO 120
7028.   C**** Check temperatures over Glacial Ice
7029.         TGI1 = (HGI(I,J,1)/XGI1/MGI(I,J) + ELHM) / SHCI
7030.         TGI2 = (HGI(I,J,2)/XGI2/MGI(I,J) + ELHM) / SHCI
7031.         TGI3 = (HGI(I,J,3)/XGI3/ACE2GI   + ELHM) / SHCI
7032.         TGI4 = (HGI(I,J,4)/XGI4/ACE2GI   + ELHM) / SHCI
7033.         IF(TGI1.le.1.D-12 .and. TGI2.le.1.D-12 .and. TGI1.ge.-90. .and.
7034.        *   TGI3.le.1.D-12 .and. TGI4.le.1.D-12)  GO TO 120
7035.         WRITE (6,*) 'After ',SUBRN,': I,J,MGI,TGI=',I,J,MGI(I,J),
7036.        *      TGI1,TGI2,TGI3,TGI4
7037.         STOP 'CHECKT'
7038.     120 IF(FGRND(I,J).le.0.)  GO TO 150
7039.   C**** Check temperature of first layer over Ground
7040.         IF(BLDATA(I,J,1).ge.-80. .and. BLDATA(I,J,1).le.75.)  GO TO 150
7041.         WRITE (6,*) 'After ',SUBRN,': I,J,TG1=',I,J,BLDATA(I,J,1)
7042.         STOP 'CHECKT'
7043.   C**** Check wetness over ground
7044.   C 130 IF(BLDATA(I,J,2).ge.0. .and. BLDATA(I,J,2).le.1.)  GO TO 150
7045.   C     WRITE (6,*) 'After ',SUBRN,': I,J,L,WEARTH=',I,J,BLDATA(I,J,2)
7046.   C     STOP 'CHECKT'
7047.     150 continue
7048.         IF(.not.QOCEAN)  RETURN
7049.   C****
7050.   C**** Check ocean values for reasonable values
7051.   C****
7052.         DO 260 J=2,JM
7053.         IMAX=IM
7054.         IF(J.eq.JM)  IMAX=1
7055.         DO 260 I=1,IMAX
7056.         IF(FOCEAN(I,J).le.0.)  GO TO 260
7057.   C**** Check potential specific enthalpy of first layer over open ocean
7058.         GO1 = G0M(I,J,1)/(MO(I,J,1)*DXYP(J))
7059.         IF(GO1.ge.-10000. .and. GO1.le.200000.)  GO TO 210
7060.         WRITE (6,*) 'After ',SUBRN,': I,J,GO1=',I,J,GO1
7061.         STOP 'CHECKT'
7062.   C**** Check first layer UO velocity
7063.     210 IF(ABS(UO(I,J,1)).le.10.)  GO TO 220
7064.         WRITE (6,*) 'After ',SUBRN,': UO(1) is too large.  I,J=',I,J
7065.         WRITE (6,*) 'MO(1)=',MO(I,J,1),MO(I+1,J,1)
7066.         WRITE (6,*) 'UO(L)=',(UO(I,J,L),L=1,LMO)
7067.         WRITE (6,*) 'VO=',VO(I,J-1,1),VO(I+1,J-1,1),VO(I,J,1),VO(I+1,J,1)
7068.         WRITE (6,*) 'G0(1)=',G0M(I,J,1)*BYDXYP(J)/MO(I,J,1),
7069.        *                     G0M(I+1,J,1)*BYDXYP(J)/MO(I+1,J,1)
7070.         WRITE (6,*) 'S0(1)=',S0M(I,J,1)*BYDXYP(J)/MO(I,J,1),
7071.        *                     S0M(I+1,J,1)*BYDXYP(J)/MO(I+1,J,1)
7072.         WRITE (6,*) 'RSI=',RSI(I,J),RSI(I+1,J)
7073.         WRITE (6,*) 'MSI=',MSI(I,J,1),MSI(I+1,J,1),MSI(I,J,2),MSI(I+1,J,2)
7074.         WRITE (6,*) 'TSI1=',(HSI(I,J,1)/MSI(I,J,1)/XSI1 + ELHM) / SHCI,
7075.        *                    (HSI(I+1,J,1)/MSI(I+1,J,1)/XSI1 + ELHM) / SHCI
7076.         WRITE (6,*) 'PSI=',PSI(I,J),PSI(I+1,J)
7077.         WRITE (6,*) 'USI=',USI(I,J)
7078.         WRITE (6,*) 'VSI=',VSI(I,J-1),VSI(I+1,J-1),VSI(I,J),VSI(I+1,J)
7079.         STOP 'CHECKT'
7080.   C**** Check first layer VO velocity
7081.     220 IF(ABS(VO(I,J,1)).le.8.)  GO TO 230
7082.         WRITE (6,*) 'After ',SUBRN,': VO(1) is too large.  I,J=',I,J
7083.         WRITE (6,*) 'MO(1)=',MO(I,J,1),MO(I,J+1,1)
7084.         WRITE (6,*) 'UO=',UO(I-1,J,1),UO(I,J,1),UO(I-1,J+1,1),UO(I,J+1,1)
7085.         WRITE (6,*) 'VO(L)=',(VO(I,J,L),L=1,5)
7086.         WRITE (6,*) 'G0(1)=',G0M(I,J,1)*BYDXYP(J)/MO(I,J,1),
7087.        *                     G0M(I,J+1,1)*BYDXYP(J)/MO(I,J+1,1)
7088.         WRITE (6,*) 'S0(1)=',S0M(I,J,1)*BYDXYP(J)/MO(I,J,1),
7089.        *                     S0M(I,J+1,1)*BYDXYP(J)/MO(I,J+1,1)
7090.         WRITE (6,*) 'RSI=',RSI(I,J),RSI(I,J+1)
7091.         WRITE (6,*) 'MSI=',MSI(I,J,1),MSI(I,J+1,1),MSI(I,J,2),MSI(I,J+1,2)
7092.         WRITE (6,*) 'TSI1=',(HSI(I,J,1)/MSI(I,J,1)/XSI1 + ELHM) / SHCI,
7093.        *                    (HSI(I,J+1,1)/MSI(I,J+1,1)/XSI1 + ELHM) / SHCI
7094.         WRITE (6,*) 'PSI=',PSI(I,J),PSI(I,J+1)
7095.         WRITE (6,*) 'USI=',USI(I-1,J),USI(I,J),USI(I-1,J+1),VSI(I,J+1)
7096.         WRITE (6,*) 'VSI=',VSI(I,J)
7097.         STOP 'CHECKT'
7098.   C**** Check first layer ocean mass
7099.     230 IF(MO(I,J,1).ge.2000.)  GO TO 240
7100.         WRITE (6,*) 'After ',SUBRN,': I,J,MO,MSI2=',I,J,MO(I,J,1),
7101.        *  MSI(I,J,2)
7102.         STOP 'CHECKT'
7103.   C**** Check ocean salinity in each eighth box for the first layer
7104.     240 SALIM = .043
7105.         IF(I.eq.47 .and. J.eq.30)  SALIM = .060
7106.         IF(.5*ABS(SXM(I,J,1))+.5*ABS(SYM(I,J,1))+RRT3*ABS(SZM(I,J,1))
7107.        *      .lt.S0M(I,J,1) .and.
7108.        *  (.5*ABS(SXM(I,J,1))+.5*ABS(SYM(I,J,1))+RRT3*ABS(SZM(I,J,1))+
7109.        *          S0M(I,J,1))/(MO(I,J,1)*DXYP(J)).lt.SALIM)  GO TO 250
7110.         WRITE (6,*) 'After ',SUBRN,': I,J,S0,SX,SY,SZ=',I,J,
7111.        *  S0M(I,J,1)/(MO(I,J,1)*DXYP(J)),SXM(I,J,1)/(MO(I,J,1)*DXYP(J)),
7112.        *  SYM(I,J,1)/(MO(I,J,1)*DXYP(J)),SZM(I,J,1)/(MO(I,J,1)*DXYP(J))
7112.1        WRITE (6,*) ' S0(L)=',(S0M(I,J,L)/(MO(I,J,L)*DXYP(J)),L=1,LMO)
7113.         STOP 'CHECKT'
7114.   C**** Check horizontal sea ice ratio to ocean
7115.     250 IF((RSI(I,J)-ABS(RSIX(I,J)).ge.  -1.d-12) .and.
7116.        *   (RSI(I,J)+ABS(RSIX(I,J)).le.1.+1.d-12) .and.
7117.        *   (RSI(I,J)-ABS(RSIY(I,J)).ge.  -1.d-12) .and.
7118.        *   (RSI(I,J)+ABS(RSIY(I,J)).le.1.+1.d-12))  GO TO 260
7119.         WRITE (6,*) 'After ',SUBRN,': I,J,RSI,RSIX,RSIY=',I,J,RSI(I,J),
7120.        *  RSIX(I,J),RSIY(I,J)
7121.         WRITE (6,*) 'USI=',USI(I-1,J),USI(I,J)
7122.         WRITE (6,*) 'VSI=',VSI(I,J-1),VSI(I,J)
7123.         WRITE (6,*) 'UA(1)=',UA(I-1,J,1),UA(I,J,1)
7124.         WRITE (6,*) 'VA(1)=',VA(I,J-1,1),VA(I,J,1)
7125.         WRITE (6,*) 'UO(1)=',UO(I-1,J,1),UO(I,J,1)
7126.         WRITE (6,*) 'VO(1)=',VO(I,J-1,1),VO(I,J,1)
7127.         STOP 'CHECKT'
7128.     260 continue
7129.         RETURN
7130.         END
8000.   
8001.         SUBROUTINE TIMER (MNOW,MSUM)
8002.   C****
8003.   C**** Output: MNOW (.01 s) = current CPU time
8004.   C****         MSUM (.01 s) = MSUM + time since last call to TIMER
8005.   C****
8006.         DATA MLAST /0/
8007.         MNOW  = MCLOCK ()
8008.         MSUM  = MNOW - MLAST + MSUM
8009.         MLAST = MNOW
8010.         RETURN
8011.   C****
8012.   C****
8013.   C**** TIMER0 and TIMER1 are used to time an individual subroutine and
8014.   C**** to attribute its cost to a short process (e.g. diagnostics) while
8015.   C**** a longer process is being timed both before and after the short
8016.   C**** subroutine was called.
8017.   C****
8018.         ENTRY TIMER0
8019.         MNOW0 = MCLOCK ()
8020.         RETURN
8021.   C****
8022.   C****
8023.         ENTRY TIMER1 (MSUM)
8024.         MNOW1 = MCLOCK ()
8025.         MSUM  = MNOW1 - MNOW0 + MSUM
8026.         MLAST = MNOW1 - MNOW0 + MLAST
8027.         RETURN
8028.         END
9000.   
9001.         SUBROUTINE UNUSED
9002.   C****
9003.   C**** Ocean subroutines that are not loaded
9004.   C****
9030.   C****
9031.   C**** Sea ice subroutines that are not needed
9032.   C****
9042.   C****
9043.   C**** Atmospheric climatology is not needed
9044.   C****
9045.         ENTRY ACLIM
9046.         RETURN
9047.         END