/[MITgcm]/MITgcm_contrib/jscott/pkg_atm2d/init_atm2d.F
ViewVC logotype

Contents of /MITgcm_contrib/jscott/pkg_atm2d/init_atm2d.F

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


Revision 1.1 - (show annotations) (download)
Fri Aug 11 18:55:50 2006 UTC (19 years, 4 months ago) by jscott
Branch: MAIN
new 2d atm package

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.COM"
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 INTEGER ib, ibj1, ibj2
46 INTEGER j_atm, mn
47 _RL end1, end2
48 _RL totrun_b(1:sNy)
49 _RL a1,a2
50 _RS atm_dyG(1:jm0)
51 DATA atm_dyG/2.0,44*4.0,2.0/
52
53 dtatmo = dtatm * 3600.
54 dtocno = dtocn * 3600.
55 dtcouplo= dtcouple * 3600.
56
57 C override data.ice seaice time step parms
58 C these will need to change if coupling procedure changed
59 thSice_deltaT = dtcouplo
60 thsIce_dtTemp = dtatmo
61 ocean_deltaT = dtcouplo
62
63 CJRS This next check - only kill it if not MPI?
64 IF (dtocno.NE.dTtracerLev(1)) THEN
65 PRINT *,'Ocean tracer timestep differs between coupler '
66 PRINT *,'and the ocean data file'
67 STOP
68 ENDIF
69
70 c Assuming the atmospheric grid array not passed, do this:
71 atm_yG(1)=-90.0
72 DO j_atm=2,jm0
73 atm_yG(j_atm)=atm_yG(j_atm-1)+atm_dyG(j_atm-1)
74 atm_yC(j_atm-1)=(atm_yG(j_atm-1)+atm_yG(j_atm))/2.0
75 ENDDO
76 atm_yC(jm0)=atm_yG(jm0)+atm_dyG(jm0)/2.0
77
78 c end atmos grid initialization
79
80 atm_oc_ind(1)=2
81 atm_oc_wgt(1)=1. _d 0
82 atm_oc_ind(sNy)=jm0-1
83 atm_oc_wgt(sNy)=1. _d 0
84
85 DO j=2, sNy-1
86
87 DO jj=2,jm0-1
88 IF ((yG(1,j,1,1).ge.atm_yG(jj)).AND.
89 & (yG(1,j,1,1).lt.atm_yG(jj+1))) j_atm=jj
90 ENDDO
91
92 atm_oc_ind(j)=j_atm
93 IF ( yG(1,j+1,1,1) .gt. atm_yG(j_atm+1) ) THEN
94 end1= sin(yG(1,j,1,1) *deg2rad)
95 end2= sin(yG(1,j+1,1,1) *deg2rad)
96 atm_oc_wgt(j)=(sin(atm_yG(j_atm+1) *deg2rad)-end1)/
97 & (end2-end1)
98 ELSE
99 atm_oc_wgt(j)=1. _d 0
100 ENDIF
101 ENDDO
102
103 c
104 c find land fraction
105 c
106 DO j_atm=1,jm0
107 cflan(j_atm)=0. _d 0
108 ocnArea(j_atm)=0. _d 0
109 ENDDO
110
111 DO j=1,sNy
112 DO i=1,sNx
113 IF (maskC(i,j,1,1,1).EQ.1.) THEN
114 ocnArea(atm_oc_ind(j))=ocnArea(atm_oc_ind(j)) +
115 & rA(i,j,1,1)*atm_oc_wgt(j)
116 IF (atm_oc_wgt(j).lt.1.d0) THEN
117 ocnArea(atm_oc_ind(j)+1)=ocnArea(atm_oc_ind(j)+1) +
118 & rA(i,j,1,1)*(1.d0-atm_oc_wgt(j))
119 ENDIF
120 ENDIF
121 ENDDO
122 ENDDO
123
124 DO j_atm=3,jm0-2
125 cflan(j_atm)=1. _d 0 - ocnArea(j_atm) /
126 & (2. _d 0 * PI * 6.37 _d 6 * 6.37 _d 6 *
127 & (sin(atm_yG(j_atm+1)*deg2rad) - sin(atm_yG(j_atm)*deg2rad)))
128 if (cflan(j_atm).LT.1. _d -14) cflan(j_atm)=0. _d 0
129 ENDDO
130
131 C deal with the combined atmos grid end cells...
132 cflan(2)= 1. _d 0 - ocnArea(2) /
133 & (2. _d 0*PI*6.37 _d 6*6.37 _d 6*
134 & (sin(atm_yG(3)*deg2rad)+1. _d 0))
135 IF (cflan(2).LT.1. _d -14) cflan(2)=0. _d 0
136 cflan(1)=cflan(2)
137 cflan(jm0-1)= 1.d0 - ocnArea(jm0-1) /
138 & (2. _d 0*PI*6.37 _d 6*6.37 _d 6*
139 & (1. _d 0-sin(atm_yG(jm0-1)*deg2rad)))
140 IF (cflan(jm0-1).LT.1. _d -14) cflan(jm0-1)=0. _d 0
141 cflan(jm0)=cflan(jm0-1)
142
143 PRINT *,'Land fractions on atmospheric grid: '
144 PRINT *, cflan
145 PRINT *,'Lookup grid index, weights:'
146 PRINT *, atm_oc_ind,atm_oc_wgt
147
148 c
149 c read in mean 1D atmos wind files -- store in memory
150 c
151 DO j_atm=1,jm0
152 DO mn=1,nForcingPer
153 atau(j_atm,mn)=0. _d 0
154 atav(j_atm,mn)=0. _d 0
155 awind(j_atm,mn)=0. _d 0
156 ENDDO
157 ENDDO
158
159 IF ( atmosTauuFile .NE. ' ' ) THEN
160 OPEN(UNIT=97, FILE=atmosTauuFile,STATUS='old',
161 & ACCESS='direct', RECL=8*jm0*nForcingPer,
162 & FORM='unformatted')
163 READ(97,REC=1), atau
164 CLOSE(97)
165 ENDIF
166
167 IF ( atmosTauvFile .NE. ' ' ) THEN
168 OPEN(UNIT=98, FILE=atmosTauvFile, STATUS='old',
169 & ACCESS='direct', RECL=8*jm0*nForcingPer,
170 & FORM='unformatted')
171 READ(98, REC=1), atav
172 CLOSE(98)
173 ENDIF
174
175 IF ( atmosWindFile .NE. ' ' ) THEN
176 OPEN(UNIT=99, FILE=atmosWindFile, STATUS='old',
177 & ACCESS='direct', RECL=8*jm0*nForcingPer,
178 & FORM='unformatted')
179 READ(99, REC=1), awind
180 CLOSE(99)
181 ENDIF
182
183 C The polar data point values for winds are effectively N/A given the
184 C pole issue... although they are read in here, they are never used in
185 C any calculations, as the polar ocean points access the data at atmos
186 C 2 and jm0-1 points.
187
188
189 c read in runoff data
190 c to put runoff into specific grid cells
191 c
192 IF ( runoffMapFile .NE. ' ' ) THEN
193 CALL READ_FLD_XY_RL( runoffMapFile, ' ',
194 & runoffVal, 0, myThid )
195 ELSE
196 DO j=1,sNy
197 DO i=1,sNx
198 if ( (maskC(i,j,1,1,1).EQ.1.) .AND.
199 & ( (maskC(i-1,j,1,1,1).EQ.0.).OR.
200 & (maskC(i+1,j,1,1,1).EQ.0.).OR.
201 & (maskC(i,j-1,1,1,1).EQ.0.).OR.
202 & (maskC(i,j+1,1,1,1).EQ.0.).OR.
203 & (maskC(i+1,j+1,1,1,1).EQ.0.).OR.
204 & (maskC(i-1,j-1,1,1,1).EQ.0.).OR.
205 & (maskC(i+1,j-1,1,1,1).EQ.0.).OR.
206 & (maskC(i-1,j+1,1,1,1).EQ.0.) ) ) THEN
207 runoffVal(i,j)=1. _d 0
208 ENDIF
209 ENDDO
210 ENDDO
211 ENDIF
212
213 DO ib=1,numBands
214 ibj1=1
215 if (ib.GT.1) ibj1= rband(ib-1)+1
216 ibj2=sNy
217 if (ib.LT.numBands) ibj2= rband(ib)
218 totrun_b(ib)=0.d0
219 DO j=ibj1,ibj2
220 DO i=1,sNx
221 totrun_b(ib)=totrun_b(ib)+runoffVal(i,j)*maskC(i,j,1,1,1)
222 ENDDO
223 ENDDO
224 DO j=ibj1,ibj2
225 runIndex(j)= ib ! for lookup of rband as fn. of latitude
226 DO i=1,sNx
227 runoffVal(i,j)=runoffVal(i,j)*maskC(i,j,1,1,1)/totrun_b(ib)
228 ENDDO
229 ENDDO
230 ENDDO
231
232 CALL INIT_SUMVARS(myThid)
233
234 C Initialize 1D diagnostic variables
235 DO j_atm=1,jm0
236 DO mn=1,nForcingPer
237 sum_tauu_ta(j_atm,mn)= 0. _d 0
238 sum_tauv_ta(j_atm,mn)= 0. _d 0
239 sum_wsocean_ta(j_atm,mn)= 0. _d 0
240 sum_ps4ocean_ta(j_atm,mn)= 0. _d 0
241 ENDDO
242 ENDDO
243
244 C Initialize 2D diagnostic variables
245 DO i=1-OLx,sNx+OLx
246 DO j=1-OLy,sNy+OLy
247 DO mn=1,nForcingPer
248 qnet_atm_ta(i,j,mn)= 0. _d 0
249 evap_atm_ta(i,j,mn)= 0. _d 0
250 precip_atm_ta(i,j,mn)= 0. _d 0
251 runoff_atm_ta(i,j,mn)= 0. _d 0
252 sum_qrel_ta(i,j,mn)= 0. _d 0
253 sum_frel_ta(i,j,mn)= 0. _d 0
254 ENDDO
255 qnet_atm(i,j)= 0. _d 0
256 evap_atm(i,j)= 0. _d 0
257 precip_atm(i,j)= 0. _d 0
258 runoff_atm(i,j)= 0. _d 0
259 sum_qrel(i,j)= 0. _d 0
260 sum_frel(i,j)= 0. _d 0
261 ENDDO
262 ENDDO
263
264 C Initialize following for safety and/or cold start
265 DO i=1-OLx,sNx+OLx
266 DO j=1-OLy,sNy+OLy
267 pass_runoff(i,j)= 0. _d 0
268 pass_qnet(i,j)= 0. _d 0
269 pass_evap(i,j)= 0. _d 0
270 pass_precip(i,j)= 0. _d 0
271 pass_fu(i,j)= 0. _d 0
272 pass_fv(i,j)= 0. _d 0
273 pass_wspeed(i,j)= 0. _d 0
274 pass_solarnet(i,j)= 0. _d 0
275 pass_slp(i,j)= 0. _d 0
276 pass_siceLoad(i,j)= 0. _d 0
277 pass_pCO2(i,j)= 0. _d 0
278 pass_prcAtm(i,j)= 0. _d 0
279 sFluxFromIce(i,j)= 0. _d 0
280 ENDDO
281 ENDDO
282
283 RETURN
284 END
285

  ViewVC Help
Powered by ViewVC 1.1.22