pkg/atm2d/init_atm2d.F

#include "ctrparam.h"
#include "ATM2D_OPTIONS.h"

C     !INTERFACE:
      SUBROUTINE INIT_ATM2D(dtatm, dtocn, dtcouple, myThid )
C     *==========================================================*
C     |  INIT_1DTO2D                                             |
C     |    This initialization routine should be run after the   |
c     |    the ocean grid/pickup have been read in.               |
c     |                                                          |
c     |  Note: grid variable indices bi,bj are hard-coded 1,1    |
c     |  This should work if coupler or atmos/coupler on one     |
c     |  machine.                                                |
c     |                                                          |
C     *==========================================================*
c
        IMPLICIT NONE

C     === Global Atmosphere Variables ===
#include "ATMSIZE.h"
#include "AGRID.h"

C     === Global Ocean Variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

C     === Global SeaIce Parameters ===
#include "THSICE_PARAMS.h"

C     === Atmos/Ocean/Seaice Interface Variables ===
#include "ATM2D_VARS.h"
 

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     dtatm, dtocn, dtcouple - Timesteps from couple.nml (hours)
C     myThid - Thread no. that called this routine.
      INTEGER dtatm, dtocn, dtcouple
      INTEGER myThid

C     LOCAL VARIABLES:
      INTEGER i,j,jj
      INTEGER ib, ibj1, ibj2  ! runoff band loop counters
      INTEGER j_atm, mn
      INTEGER dUnit
      _RL end1, end2, enda1, enda2, enda3 !used to compute grid conv areas
      _RL totrun_b(sNy) ! total file "runoff" in runoff bands
      _RL a1,a2
      _RS atm_dyG(jm0)  ! southern point/(boundary) of atmos grid
      DATA atm_dyG/2.0,44*4.0,2.0/  ! grid spacing for atmosphere

      dtatmo = dtatm * 3600.
      dtocno = dtocn * 3600.
      dtcouplo= dtcouple * 3600.

C override data.ice seaice time step parms
C these will need to change if coupling procedure changed
      thSice_deltaT = dtcouplo
      thsIce_dtTemp = dtatmo
      ocean_deltaT = dtcouplo

CJRS  This next check - only kill it if not MPI?
      IF (dtocno.NE.dTtracerLev(1)) THEN
        PRINT *,'Ocean tracer timestep differs between coupler '
        PRINT *,'and the ocean data file'
        STOP
      ENDIF

c Assuming the atmospheric grid array not passed, do this:
      atm_yG(1)=-90.0
      DO j_atm=2,jm0
        atm_yG(j_atm)=atm_yG(j_atm-1)+atm_dyG(j_atm-1)
        atm_yC(j_atm-1)=(atm_yG(j_atm-1)+atm_yG(j_atm))/2.0
      ENDDO
      atm_yC(jm0)=atm_yG(jm0)+atm_dyG(jm0)/2.0

c end atmos grid initialization

      atm_oc_ind(1)=2
      atm_oc_wgt(1)=1. _d 0
      atm_oc_frac1(1)= (sin(yG(1,2,1,1)*deg2rad) - 
     &        sin(yG(1,1,1,1)*deg2rad))/
     &        (sin(atm_yG(3)*deg2rad)-sin(atm_yG(1)*deg2rad))
      atm_oc_frac2(1)= 0. _d 0   ! assumes ocean(1) fits in atm(1)
      atm_oc_ind(sNy)=jm0-1
      atm_oc_wgt(sNy)=1. _d 0
      atm_oc_frac1(sNy)= (sin((yG(1,sNy,1,1) + 
     &      dyG(1,sNy,1,1)/6.37D6/deg2rad)*deg2rad)-
     &      sin(yG(1,sNy,1,1)*deg2rad))/
     &      (sin((atm_yG(jm0)+atm_dyG(jm0))*deg2rad)- 
     &      sin(atm_yG(jm0-1)*deg2rad))
      atm_oc_frac2(sNy)= 0. _d 0   ! assumes ocean(1) fits in atm(1)

      endwgt1 = sin(atm_yG(2)*deg2rad)            !hard-coded that the atmos
      endwgt2 = sin(atm_yG(3)*deg2rad) - endwgt1  !grid is same in NH and SH
      endwgt1 = endwgt1 + 1. _d 0                 !and goes 90S to 90N
      rsumwgt = 1. _d 0/(endwgt1 + endwgt2)

      atm_yG(2)=atm_yG(1) ! grid now combined atm end points
      atm_yG(jm0)=90. _d 0

      DO j=2, sNy-1

        DO jj=2,jm0-1
          IF  ((yG(1,j,1,1).GE.atm_yG(jj)).AND.
     &         (yG(1,j,1,1).LT.atm_yG(jj+1))) j_atm=jj
        ENDDO  
        
        atm_oc_ind(j)=j_atm
        end1= sin(yG(1,j,1,1) *deg2rad)
        end2= sin(yG(1,j+1,1,1) *deg2rad)
        enda1 = sin(atm_yG(j_atm) *deg2rad)
        enda2 = sin(atm_yG(j_atm+1) *deg2rad)
        IF ( yG(1,j+1,1,1) .GT. atm_yG(j_atm+1) ) THEN
           enda3 = sin(atm_yG(j_atm+2) *deg2rad)
           atm_oc_wgt(j)=(enda2-end1)/ (end2-end1)
           atm_oc_frac1(j)= (enda2-end1) / (enda2 - enda1)
           atm_oc_frac2(j)= (end2 - enda2) / (enda3 - enda2)
        ELSE
          atm_oc_wgt(j)=1. _d 0
          atm_oc_frac1(j)= (end2-end1)/ (enda2-enda1)
          atm_oc_frac2(j)=0. _d 0
        ENDIF
      ENDDO
c
c find land fraction
c
      DO j_atm=1,jm0
        cflan(j_atm)=0. _d 0
        ocnArea(j_atm)=0. _d 0
      ENDDO

      DO j=1,sNy
        DO i=1,sNx
          IF (maskC(i,j,1,1,1).EQ.1.) THEN
            ocnArea(atm_oc_ind(j))=ocnArea(atm_oc_ind(j)) +
     &                           rA(i,j,1,1)*atm_oc_wgt(j)
            IF (atm_oc_wgt(j).LT.1.d0) THEN
              ocnArea(atm_oc_ind(j)+1)=ocnArea(atm_oc_ind(j)+1) +
     &                           rA(i,j,1,1)*(1.d0-atm_oc_wgt(j))
            ENDIF
          ENDIF
        ENDDO
      ENDDO

      DO j_atm=3,jm0-2
        cflan(j_atm)=1. _d 0 - ocnArea(j_atm) / 
     &              (2. _d 0 * PI * 6.37 _d 6 * 6.37 _d 6 * 
     &    (sin(atm_yG(j_atm+1)*deg2rad) - sin(atm_yG(j_atm)*deg2rad)))
        if (cflan(j_atm).LT.1. _d -14) cflan(j_atm)=0. _d 0
      ENDDO

C     deal with the combined atmos grid end cells...
      cflan(2)= 1. _d 0 - ocnArea(2) / 
     &         (2. _d 0*PI*6.37 _d 6*6.37 _d 6*
     &         (sin(atm_yG(3)*deg2rad)+1. _d 0))
      IF (cflan(2).LT.1. _d -14) cflan(2)=0. _d 0
      cflan(1)=cflan(2)
      cflan(jm0-1)= 1.d0 - ocnArea(jm0-1) / 
     &             (2. _d 0*PI*6.37 _d 6*6.37 _d 6*
     &             (1. _d 0-sin(atm_yG(jm0-1)*deg2rad)))
      IF (cflan(jm0-1).LT.1. _d -14) cflan(jm0-1)=0. _d 0
      cflan(jm0)=cflan(jm0-1)

      PRINT *,'Land fractions on atmospheric grid: ' 
      PRINT *, cflan
      PRINT *,'Lookup grid index, weights:'
      PRINT *, atm_oc_ind,atm_oc_wgt
C      PRINT *,'Lookup fraction 1 of atmos grid:'
C      PRINT *, atm_oc_frac1
C      PRINT *,'Lookup fraction 2 of atmos grid:'
C      PRINT *, atm_oc_frac2

c
c read in mean 1D atmos wind files -- store in memory
c
      DO j_atm=1,jm0
        DO mn=1,nForcingPer
          atau(j_atm,mn)=0. _d 0
          atav(j_atm,mn)=0. _d 0
          awind(j_atm,mn)=0. _d 0
        ENDDO
      ENDDO

      CALL MDSFINDUNIT( dUnit, myThid )

      IF ( atmosTauuFile .NE. ' '  ) THEN
         OPEN(dUnit, FILE=atmosTauuFile,STATUS='old',
     &        ACCESS='direct', RECL=8*jm0*nForcingPer, 
     &        FORM='unformatted')
         READ(dUnit,REC=1), atau
         CLOSE(dUnit)
      ENDIF

      IF ( atmosTauvFile .NE. ' '  ) THEN
         OPEN(dUnit, FILE=atmosTauvFile, STATUS='old',
     &        ACCESS='direct', RECL=8*jm0*nForcingPer, 
     &        FORM='unformatted')
         READ(dUnit, REC=1), atav
         CLOSE(dUnit)
      ENDIF

      IF ( atmosWindFile .NE. ' '  ) THEN
         OPEN(dUnit, FILE=atmosWindFile, STATUS='old',
     &        ACCESS='direct', RECL=8*jm0*nForcingPer, 
     &        FORM='unformatted')
         READ(dUnit, REC=1), awind
         CLOSE(dUnit)
      ENDIF

C The polar data point values for winds are effectively N/A given the
C pole issue... although they are read in here, they are never used in
C any calculations, as the polar ocean points access the data at atmos
C 2 and jm0-1 points.


c read in runoff data
c to put runoff into specific grid cells
c
      IF ( runoffMapFile .NE. ' ' ) THEN
        CALL READ_FLD_XY_RL( runoffMapFile, ' ',
     &                      runoffVal, 0, myThid )
      ELSE
        DO j=1,sNy
          DO i=1,sNx
            if ( (maskC(i,j,1,1,1).EQ.1.) .AND.
     &          ( (maskC(i-1,j,1,1,1).EQ.0.).OR.
     &            (maskC(i+1,j,1,1,1).EQ.0.).OR.
     &            (maskC(i,j-1,1,1,1).EQ.0.).OR.
     &            (maskC(i,j+1,1,1,1).EQ.0.).OR.
     &            (maskC(i+1,j+1,1,1,1).EQ.0.).OR.
     &            (maskC(i-1,j-1,1,1,1).EQ.0.).OR.
     &            (maskC(i+1,j-1,1,1,1).EQ.0.).OR.
     &            (maskC(i-1,j+1,1,1,1).EQ.0.) ) ) THEN
              runoffVal(i,j)=1. _d 0
            ENDIF
          ENDDO
        ENDDO
      ENDIF

      DO ib=1,numBands
        ibj1=1
        if (ib.GT.1) ibj1=  rband(ib-1)+1
        ibj2=sNy
        if (ib.LT.numBands) ibj2= rband(ib)
        totrun_b(ib)=0.d0
        DO j=ibj1,ibj2
          DO i=1,sNx
            totrun_b(ib)=totrun_b(ib)+runoffVal(i,j)*maskC(i,j,1,1,1)
          ENDDO
        ENDDO
        DO j=ibj1,ibj2
          runIndex(j)= ib     ! for lookup of rband as fn. of latitude
          DO i=1,sNx
            runoffVal(i,j)=runoffVal(i,j)*maskC(i,j,1,1,1)/totrun_b(ib)
          ENDDO
        ENDDO
      ENDDO

      CALL INIT_SUMVARS(myThid)
 
C     Initialize 1D diagnostic variables 
      DO j_atm=1,jm0
        DO mn=1,nForcingPer
          sum_tauu_ta(j_atm,mn)= 0. _d 0
          sum_tauv_ta(j_atm,mn)= 0. _d 0
          sum_wsocean_ta(j_atm,mn)= 0. _d 0
          sum_ps4ocean_ta(j_atm,mn)= 0. _d 0
        ENDDO
      ENDDO

C     Initialize 2D diagnostic variables 
      DO i=1-OLx,sNx+OLx
        DO j=1-OLy,sNy+OLy
          DO mn=1,nForcingPer
            qnet_atm_ta(i,j,mn)= 0. _d 0
            evap_atm_ta(i,j,mn)= 0. _d 0
            precip_atm_ta(i,j,mn)= 0. _d 0
            runoff_atm_ta(i,j,mn)= 0. _d 0
            sum_qrel_ta(i,j,mn)= 0. _d 0
            sum_frel_ta(i,j,mn)= 0. _d 0
            sum_iceMask_ta(i,j,mn)= 0. _d 0
            sum_iceHeight_ta(i,j,mn)= 0. _d 0
            sum_iceTime_ta(i,j,mn)= 0. _d 0
            sum_oceMxLT_ta(i,j,mn)= 0. _d 0
            sum_oceMxLS_ta(i,j,mn)= 0. _d 0
          ENDDO
          qnet_atm(i,j)= 0. _d 0
          evap_atm(i,j)= 0. _d 0
          precip_atm(i,j)= 0. _d 0
          runoff_atm(i,j)= 0. _d 0
          sum_qrel(i,j)= 0. _d 0
          sum_frel(i,j)= 0. _d 0
          sum_iceMask(i,j)= 0. _d 0
          sum_iceHeight(i,j)= 0. _d 0
          sum_iceTime(i,j)= 0. _d 0
          sum_oceMxLT(i,j)= 0. _d 0
          sum_oceMxLS(i,j)= 0. _d 0
        ENDDO
      ENDDO

C     Initialize following for safety and/or cold start
      DO i=1-OLx,sNx+OLx
        DO j=1-OLy,sNy+OLy
          pass_runoff(i,j)= 0. _d 0
          pass_qnet(i,j)= 0. _d 0
          pass_evap(i,j)= 0. _d 0
          pass_precip(i,j)= 0. _d 0
          pass_fu(i,j)= 0. _d 0
          pass_fv(i,j)= 0. _d 0
          pass_wspeed(i,j)= 0. _d 0
          pass_solarnet(i,j)= 0. _d 0
          pass_slp(i,j)= 0. _d 0
          pass_siceLoad(i,j)= 0. _d 0
          pass_pCO2(i,j)= 0. _d 0
          pass_prcAtm(i,j)= 0. _d 0
          sFluxFromIce(i,j)= 0. _d 0
        ENDDO
      ENDDO

      RETURN
      END

1	jscott	1.1	#include "ctrparam.h"
2			#include "ATM2D_OPTIONS.h"
3
4			C !INTERFACE:
5			SUBROUTINE INIT_ATM2D(dtatm, dtocn, dtcouple, myThid )
6			C ==========================================================
7			C \| INIT_1DTO2D \|
8			C \| This initialization routine should be run after the \|
9			c \| the ocean grid/pickup have been read in. \|
10			c \| \|
11			c \| Note: grid variable indices bi,bj are hard-coded 1,1 \|
12			c \| This should work if coupler or atmos/coupler on one \|
13			c \| machine. \|
14			c \| \|
15			C ==========================================================
16			c
17			IMPLICIT NONE
18
19			C === Global Atmosphere Variables ===
20			#include "ATMSIZE.h"
21			#include "AGRID.h"
22
23			C === Global Ocean Variables ===
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "GRID.h"
28
29			C === Global SeaIce Parameters ===
30			#include "THSICE_PARAMS.h"
31
32			C === Atmos/Ocean/Seaice Interface Variables ===
33			#include "ATM2D_VARS.h"
34
35
36			C !INPUT/OUTPUT PARAMETERS:
37			C === Routine arguments ===
38			C dtatm, dtocn, dtcouple - Timesteps from couple.nml (hours)
39			C myThid - Thread no. that called this routine.
40			INTEGER dtatm, dtocn, dtcouple
41			INTEGER myThid
42
43			C LOCAL VARIABLES:
44			INTEGER i,j,jj
45	jscott	1.2	INTEGER ib, ibj1, ibj2 ! runoff band loop counters
46	jscott	1.1	INTEGER j_atm, mn
47			INTEGER dUnit
48	jscott	1.2	_RL end1, end2, enda1, enda2, enda3 !used to compute grid conv areas
49			_RL totrun_b(sNy) ! total file "runoff" in runoff bands
50	jscott	1.1	_RL a1,a2
51	jscott	1.2	_RS atm_dyG(jm0) ! southern point/(boundary) of atmos grid
52			DATA atm_dyG/2.0,44*4.0,2.0/ ! grid spacing for atmosphere
53	jscott	1.1
54			dtatmo = dtatm * 3600.
55			dtocno = dtocn * 3600.
56			dtcouplo= dtcouple * 3600.
57
58			C override data.ice seaice time step parms
59			C these will need to change if coupling procedure changed
60			thSice_deltaT = dtcouplo
61			thsIce_dtTemp = dtatmo
62			ocean_deltaT = dtcouplo
63
64			CJRS This next check - only kill it if not MPI?
65			IF (dtocno.NE.dTtracerLev(1)) THEN
66			PRINT *,'Ocean tracer timestep differs between coupler '
67			PRINT *,'and the ocean data file'
68			STOP
69			ENDIF
70
71			c Assuming the atmospheric grid array not passed, do this:
72			atm_yG(1)=-90.0
73			DO j_atm=2,jm0
74			atm_yG(j_atm)=atm_yG(j_atm-1)+atm_dyG(j_atm-1)
75			atm_yC(j_atm-1)=(atm_yG(j_atm-1)+atm_yG(j_atm))/2.0
76			ENDDO
77			atm_yC(jm0)=atm_yG(jm0)+atm_dyG(jm0)/2.0
78
79			c end atmos grid initialization
80
81			atm_oc_ind(1)=2
82			atm_oc_wgt(1)=1. _d 0
83	jscott	1.2	atm_oc_frac1(1)= (sin(yG(1,2,1,1)*deg2rad) -
84			& sin(yG(1,1,1,1)*deg2rad))/
85			& (sin(atm_yG(3)deg2rad)-sin(atm_yG(1)deg2rad))
86			atm_oc_frac2(1)= 0. _d 0 ! assumes ocean(1) fits in atm(1)
87	jscott	1.1	atm_oc_ind(sNy)=jm0-1
88			atm_oc_wgt(sNy)=1. _d 0
89	jscott	1.2	atm_oc_frac1(sNy)= (sin((yG(1,sNy,1,1) +
90			& dyG(1,sNy,1,1)/6.37D6/deg2rad)*deg2rad)-
91			& sin(yG(1,sNy,1,1)*deg2rad))/
92			& (sin((atm_yG(jm0)+atm_dyG(jm0))*deg2rad)-
93			& sin(atm_yG(jm0-1)*deg2rad))
94			atm_oc_frac2(sNy)= 0. _d 0 ! assumes ocean(1) fits in atm(1)
95
96			endwgt1 = sin(atm_yG(2)*deg2rad) !hard-coded that the atmos
97			endwgt2 = sin(atm_yG(3)*deg2rad) - endwgt1 !grid is same in NH and SH
98			endwgt1 = endwgt1 + 1. _d 0 !and goes 90S to 90N
99			rsumwgt = 1. _d 0/(endwgt1 + endwgt2)
100
101			atm_yG(2)=atm_yG(1) ! grid now combined atm end points
102			atm_yG(jm0)=90. _d 0
103	jscott	1.1
104			DO j=2, sNy-1
105
106			DO jj=2,jm0-1
107	jscott	1.2	IF ((yG(1,j,1,1).GE.atm_yG(jj)).AND.
108			& (yG(1,j,1,1).LT.atm_yG(jj+1))) j_atm=jj
109	jscott	1.1	ENDDO
110
111			atm_oc_ind(j)=j_atm
112	jscott	1.2	end1= sin(yG(1,j,1,1) *deg2rad)
113			end2= sin(yG(1,j+1,1,1) *deg2rad)
114			enda1 = sin(atm_yG(j_atm) *deg2rad)
115			enda2 = sin(atm_yG(j_atm+1) *deg2rad)
116			IF ( yG(1,j+1,1,1) .GT. atm_yG(j_atm+1) ) THEN
117			enda3 = sin(atm_yG(j_atm+2) *deg2rad)
118			atm_oc_wgt(j)=(enda2-end1)/ (end2-end1)
119			atm_oc_frac1(j)= (enda2-end1) / (enda2 - enda1)
120			atm_oc_frac2(j)= (end2 - enda2) / (enda3 - enda2)
121	jscott	1.1	ELSE
122			atm_oc_wgt(j)=1. _d 0
123	jscott	1.2	atm_oc_frac1(j)= (end2-end1)/ (enda2-enda1)
124			atm_oc_frac2(j)=0. _d 0
125	jscott	1.1	ENDIF
126			ENDDO
127			c
128			c find land fraction
129			c
130			DO j_atm=1,jm0
131			cflan(j_atm)=0. _d 0
132			ocnArea(j_atm)=0. _d 0
133			ENDDO
134
135			DO j=1,sNy
136			DO i=1,sNx
137			IF (maskC(i,j,1,1,1).EQ.1.) THEN
138			ocnArea(atm_oc_ind(j))=ocnArea(atm_oc_ind(j)) +
139			& rA(i,j,1,1)*atm_oc_wgt(j)
140	jscott	1.2	IF (atm_oc_wgt(j).LT.1.d0) THEN
141	jscott	1.1	ocnArea(atm_oc_ind(j)+1)=ocnArea(atm_oc_ind(j)+1) +
142			& rA(i,j,1,1)*(1.d0-atm_oc_wgt(j))
143			ENDIF
144			ENDIF
145			ENDDO
146			ENDDO
147
148			DO j_atm=3,jm0-2
149			cflan(j_atm)=1. _d 0 - ocnArea(j_atm) /
150			& (2. _d 0 * PI * 6.37 _d 6 * 6.37 _d 6 *
151			& (sin(atm_yG(j_atm+1)deg2rad) - sin(atm_yG(j_atm)deg2rad)))
152			if (cflan(j_atm).LT.1. _d -14) cflan(j_atm)=0. _d 0
153			ENDDO
154
155			C deal with the combined atmos grid end cells...
156			cflan(2)= 1. _d 0 - ocnArea(2) /
157			& (2. _d 0PI6.37 _d 66.37 _d 6
158			& (sin(atm_yG(3)*deg2rad)+1. _d 0))
159			IF (cflan(2).LT.1. _d -14) cflan(2)=0. _d 0
160			cflan(1)=cflan(2)
161			cflan(jm0-1)= 1.d0 - ocnArea(jm0-1) /
162			& (2. _d 0PI6.37 _d 66.37 _d 6
163			& (1. _d 0-sin(atm_yG(jm0-1)*deg2rad)))
164			IF (cflan(jm0-1).LT.1. _d -14) cflan(jm0-1)=0. _d 0
165			cflan(jm0)=cflan(jm0-1)
166
167			PRINT *,'Land fractions on atmospheric grid: '
168			PRINT *, cflan
169			PRINT *,'Lookup grid index, weights:'
170			PRINT *, atm_oc_ind,atm_oc_wgt
171	jscott	1.2	C PRINT *,'Lookup fraction 1 of atmos grid:'
172			C PRINT *, atm_oc_frac1
173			C PRINT *,'Lookup fraction 2 of atmos grid:'
174			C PRINT *, atm_oc_frac2
175	jscott	1.1
176			c
177			c read in mean 1D atmos wind files -- store in memory
178			c
179			DO j_atm=1,jm0
180			DO mn=1,nForcingPer
181			atau(j_atm,mn)=0. _d 0
182			atav(j_atm,mn)=0. _d 0
183			awind(j_atm,mn)=0. _d 0
184			ENDDO
185			ENDDO
186
187			CALL MDSFINDUNIT( dUnit, myThid )
188
189			IF ( atmosTauuFile .NE. ' ' ) THEN
190			OPEN(dUnit, FILE=atmosTauuFile,STATUS='old',
191			& ACCESS='direct', RECL=8jm0nForcingPer,
192			& FORM='unformatted')
193			READ(dUnit,REC=1), atau
194			CLOSE(dUnit)
195			ENDIF
196
197			IF ( atmosTauvFile .NE. ' ' ) THEN
198			OPEN(dUnit, FILE=atmosTauvFile, STATUS='old',
199			& ACCESS='direct', RECL=8jm0nForcingPer,
200			& FORM='unformatted')
201			READ(dUnit, REC=1), atav
202			CLOSE(dUnit)
203			ENDIF
204
205			IF ( atmosWindFile .NE. ' ' ) THEN
206			OPEN(dUnit, FILE=atmosWindFile, STATUS='old',
207			& ACCESS='direct', RECL=8jm0nForcingPer,
208			& FORM='unformatted')
209			READ(dUnit, REC=1), awind
210			CLOSE(dUnit)
211			ENDIF
212
213			C The polar data point values for winds are effectively N/A given the
214			C pole issue... although they are read in here, they are never used in
215			C any calculations, as the polar ocean points access the data at atmos
216			C 2 and jm0-1 points.
217
218
219			c read in runoff data
220			c to put runoff into specific grid cells
221			c
222			IF ( runoffMapFile .NE. ' ' ) THEN
223			CALL READ_FLD_XY_RL( runoffMapFile, ' ',
224			& runoffVal, 0, myThid )
225			ELSE
226			DO j=1,sNy
227			DO i=1,sNx
228			if ( (maskC(i,j,1,1,1).EQ.1.) .AND.
229			& ( (maskC(i-1,j,1,1,1).EQ.0.).OR.
230			& (maskC(i+1,j,1,1,1).EQ.0.).OR.
231			& (maskC(i,j-1,1,1,1).EQ.0.).OR.
232			& (maskC(i,j+1,1,1,1).EQ.0.).OR.
233			& (maskC(i+1,j+1,1,1,1).EQ.0.).OR.
234			& (maskC(i-1,j-1,1,1,1).EQ.0.).OR.
235			& (maskC(i+1,j-1,1,1,1).EQ.0.).OR.
236			& (maskC(i-1,j+1,1,1,1).EQ.0.) ) ) THEN
237			runoffVal(i,j)=1. _d 0
238			ENDIF
239			ENDDO
240			ENDDO
241			ENDIF
242
243			DO ib=1,numBands
244			ibj1=1
245			if (ib.GT.1) ibj1= rband(ib-1)+1
246			ibj2=sNy
247			if (ib.LT.numBands) ibj2= rband(ib)
248			totrun_b(ib)=0.d0
249			DO j=ibj1,ibj2
250			DO i=1,sNx
251			totrun_b(ib)=totrun_b(ib)+runoffVal(i,j)*maskC(i,j,1,1,1)
252			ENDDO
253			ENDDO
254			DO j=ibj1,ibj2
255			runIndex(j)= ib ! for lookup of rband as fn. of latitude
256			DO i=1,sNx
257			runoffVal(i,j)=runoffVal(i,j)*maskC(i,j,1,1,1)/totrun_b(ib)
258			ENDDO
259			ENDDO
260			ENDDO
261
262			CALL INIT_SUMVARS(myThid)
263
264			C Initialize 1D diagnostic variables
265			DO j_atm=1,jm0
266			DO mn=1,nForcingPer
267			sum_tauu_ta(j_atm,mn)= 0. _d 0
268			sum_tauv_ta(j_atm,mn)= 0. _d 0
269			sum_wsocean_ta(j_atm,mn)= 0. _d 0
270			sum_ps4ocean_ta(j_atm,mn)= 0. _d 0
271			ENDDO
272			ENDDO
273
274			C Initialize 2D diagnostic variables
275			DO i=1-OLx,sNx+OLx
276			DO j=1-OLy,sNy+OLy
277			DO mn=1,nForcingPer
278			qnet_atm_ta(i,j,mn)= 0. _d 0
279			evap_atm_ta(i,j,mn)= 0. _d 0
280			precip_atm_ta(i,j,mn)= 0. _d 0
281			runoff_atm_ta(i,j,mn)= 0. _d 0
282			sum_qrel_ta(i,j,mn)= 0. _d 0
283			sum_frel_ta(i,j,mn)= 0. _d 0
284			sum_iceMask_ta(i,j,mn)= 0. _d 0
285			sum_iceHeight_ta(i,j,mn)= 0. _d 0
286			sum_iceTime_ta(i,j,mn)= 0. _d 0
287			sum_oceMxLT_ta(i,j,mn)= 0. _d 0
288			sum_oceMxLS_ta(i,j,mn)= 0. _d 0
289			ENDDO
290			qnet_atm(i,j)= 0. _d 0
291			evap_atm(i,j)= 0. _d 0
292			precip_atm(i,j)= 0. _d 0
293			runoff_atm(i,j)= 0. _d 0
294			sum_qrel(i,j)= 0. _d 0
295			sum_frel(i,j)= 0. _d 0
296			sum_iceMask(i,j)= 0. _d 0
297			sum_iceHeight(i,j)= 0. _d 0
298			sum_iceTime(i,j)= 0. _d 0
299			sum_oceMxLT(i,j)= 0. _d 0
300			sum_oceMxLS(i,j)= 0. _d 0
301			ENDDO
302			ENDDO
303
304			C Initialize following for safety and/or cold start
305			DO i=1-OLx,sNx+OLx
306			DO j=1-OLy,sNy+OLy
307			pass_runoff(i,j)= 0. _d 0
308			pass_qnet(i,j)= 0. _d 0
309			pass_evap(i,j)= 0. _d 0
310			pass_precip(i,j)= 0. _d 0
311			pass_fu(i,j)= 0. _d 0
312			pass_fv(i,j)= 0. _d 0
313			pass_wspeed(i,j)= 0. _d 0
314			pass_solarnet(i,j)= 0. _d 0
315			pass_slp(i,j)= 0. _d 0
316			pass_siceLoad(i,j)= 0. _d 0
317			pass_pCO2(i,j)= 0. _d 0
318			pass_prcAtm(i,j)= 0. _d 0
319			sFluxFromIce(i,j)= 0. _d 0
320			ENDDO
321			ENDDO
322
323			RETURN
324			END
325