/[MITgcm]/MITgcm/pkg/atm2d/init_atm2d.F
ViewVC logotype

Annotation of /MITgcm/pkg/atm2d/init_atm2d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.6 - (hide annotations) (download)
Mon Nov 19 22:57:31 2007 UTC (16 years, 6 months ago) by jscott
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59k, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q
Changes since 1.5: +19 -1 lines
interpolate meridional wind stress to match up with ocean C grid
CV: ----------------------------------------------------------------------

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

  ViewVC Help
Powered by ViewVC 1.1.22