/[MITgcm]/MITgcm/pkg/aim/aim_do_atmos_physics.F
ViewVC logotype

Annotation of /MITgcm/pkg/aim/aim_do_atmos_physics.F

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


Revision 1.6 - (hide annotations) (download)
Wed Aug 15 13:50:43 2001 UTC (22 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint40pre8
Changes since 1.5: +5 -6 lines
allow to run AIM with partial cell.

1 jmc 1.6 C $Header: /u/gcmpack/models/MITgcmUV/pkg/aim/aim_do_atmos_physics.F,v 1.5 2001/06/18 17:39:58 cnh Exp $
2 jmc 1.3 C $Name: $
3 adcroft 1.2
4     #include "AIM_OPTIONS.h"
5 cnh 1.4 #undef OLD_AIM_GRIG_MAPPING
6 adcroft 1.2
7 cnh 1.4 SUBROUTINE AIM_DO_ATMOS_PHYSICS( phi_hyd,
8     I bi, bj,
9     I currentTime, myThid )
10 adcroft 1.2 C /==================================================================\
11     C | S/R AIM_DO_ATMOS_PHYSICS |
12     C |==================================================================|
13     C | Interface interface between atmospheric physics package and the |
14     C | dynamical model. |
15     C | Routine calls physics pacakge after mapping model variables to |
16     C | the package grid. Package should derive and set tendency terms |
17     C | which can be included as external forcing terms in the dynamical |
18     C | tendency routines. Packages should communicate this information |
19     C | through common blocks. |
20     C \==================================================================/
21 cnh 1.4 IMPLICIT rEAL*8 (A-H,O-Z)
22 adcroft 1.2
23     C -------------- Global variables ------------------------------------
24     C Physics package
25     #include "atparam.h"
26     #include "atparam1.h"
27     INTEGER NGP
28     INTEGER NLON
29     INTEGER NLAT
30     INTEGER NLEV
31     PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
32 cnh 1.4
33 adcroft 1.2 C MITgcm
34     #include "EEPARAMS.h"
35     #include "PARAMS.h"
36     #include "DYNVARS.h"
37     #include "GRID.h"
38 cnh 1.4 #include "SURFACE.h"
39     #include "AIM_FFIELDS.h"
40    
41     C Physics package
42     #include "com_physvar.h"
43     #include "com_forcing1.h"
44     #include "Lev_def.h"
45 adcroft 1.2
46     C -------------- Routine arguments -----------------------------------
47 cnh 1.4 _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48 adcroft 1.2 _RL currentTime
49 cnh 1.4 INTEGER myThid
50     INTEGER bi, bj
51 adcroft 1.2
52     #ifdef ALLOW_AIM
53     C -------------- Local variables -------------------------------------
54     C I,J,K,I2,J2 - Loop counters
55     C tYear - Fraction into year
56     C mnthIndex - Current month
57     C prevMnthIndex - Month last time this routine was called.
58     C tmp4 - I/O buffer ( 32-bit precision )
59     C fNam - Work space for file names
60     C mnthNam - Month strings
61     C hInital - Initial height of pressure surfaces (m)
62     C pSurfs - Pressure surfaces (Pa)
63     C Katm - Atmospheric K index
64     INTEGER I
65     INTEGER I2
66     INTEGER J
67     INTEGER J2
68     INTEGER K
69     INTEGER IG0
70     INTEGER JG0
71     REAL tYear
72     INTEGER mnthIndex
73     INTEGER prevMnthIndex
74     DATA prevMnthIndex / 0 /
75     SAVE prevMnthIndex
76     LOGICAL FirstCall
77     DATA FirstCall /.TRUE./
78     SAVE FirstCall
79     LOGICAL CALLFirst
80     DATA CALLFirst /.TRUE./
81     SAVE CALLFirst
82     INTEGER nxIo
83     INTEGER nyIo
84     PARAMETER ( nxIo = 128, nyIo = 64 )
85     Real*4 tmp4(nxIo,nyIo)
86     CHARACTER*16 fNam
87     CHARACTER*3 mnthNam(12)
88     DATA mnthNam /
89     & 'jan', 'feb', 'mar', 'apr', 'may', 'jun',
90     & 'jul', 'aug', 'sep', 'oct', 'nov', 'dec' /
91     SAVE mnthNam
92     REAL hInitial(Nr)
93     REAL hInitialW(Nr)
94     DATA hInitial / 418.038,2038.54,5296.88,10090.02,17338.0/
95     SAVE hInitial
96     DATA hInitialW / 0., 1657.54, 4087.75, 8050.96,15090.4 /
97     REAL pSurfs(Nr)
98     DATA pSurfs / 950.D2,775.D2, 500.D2, 250.D2, 75.D2 /
99     SAVE pSurfs
100     REAL pSurfw(Nr)
101     DATA pSurfw /1000.D2, 900.D2, 650.D2, 350.D2, 150.D2 /
102     SAVE pSurfw
103     REAL RD
104     REAL CPAIR
105     REAL RhoG1(sNx*sNy,Nr)
106     INTEGER npasdt
107     DATA npasdt /0/
108     SAVE npasdt
109     REAL Soilqmax
110     REAL phiTotal(sNx,sNy,Nr)
111     _RL phiTCount
112     _RL phiTSum
113     _RL ans
114     real pvoltotNiv5
115     SAVE pvoltotNiv5
116     real ptotalNiv5
117     INTEGER Katm
118 cnh 1.4
119 adcroft 1.2 C
120     pGround = 1.D5
121     CPAIR = 1004
122     RD = 287
123    
124 cnh 1.4 CALL AIM_DYN2AIM( bi, bj, currentTime, myThid )
125    
126 adcroft 1.2 C Assume only one tile per proc. for now
127 cnh 1.4 IG0 = myXGlobalLo+(bi-1)*sNx
128     JG0 = myYGlobalLo+(bj-1)*sNy
129 adcroft 1.2
130     C
131     C Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
132     C Internal index mapping is linear in X and Y with a second
133     C dimension for the vertical.
134    
135     C Adjustment for heave due to mean heating/cooling
136     C ( I don't think the old formula was strictly "correct" for orography
137     C but I have implemented it as was for now. As implemented
138     C the mean heave of the bottom (K=Nr) level is calculated rather than
139     C the mean heave of the base of the atmosphere. )
140     phiTCount = 0.
141     phiTSum = 0.
142     DO K=1,Nr
143     DO J=1,sNy
144     DO I=1,sNx
145 jmc 1.3 phiTotal(I,J,K) = etaN(i,j,bi,bj)
146 adcroft 1.2 phiTCount = phiTCount + hFacC(i,j,Nr,bi,bj)
147     ENDDO
148     ENDDO
149     ENDDO
150     DO K=1,Nr
151     DO J=1,sNy
152     DO I=1,sNx
153     phiTotal(I,J,K) = phiTotal(I,J,K) +
154 cnh 1.4 & recip_rhoConst*(phi_hyd(i,j,k))
155 adcroft 1.2 ENDDO
156     ENDDO
157     ENDDO
158     DO J=1,sNy
159     DO I=1,sNx
160     phiTSum = phiTSum + phiTotal(I,J,Nr)
161     ENDDO
162     ENDDO
163     ans = phiTCount
164     C _GLOBAL_SUM_R8( phiTCount, myThid )
165     phiTcount = ans
166     ans = phiTSum
167     C _GLOBAL_SUM_R8( phiTSum, myThid )
168     phiTSum = ans
169     C ptotalniv5=phiTSum/phiTCount
170     ptotalniv5=0.
171    
172 cnh 1.4 c_jmc: Because AIM physics LSC is not applied in the stratosphere (top level),
173     c ==> move water wapor from the stratos to the surface level.
174     DO J = 1-Oly, sNy+Oly
175     DO I = 1-Olx, sNx+Olx
176     c k = k_surf(i,j,bi,bj)
177     c salt(I,J,k,bi,bj) = salt(I,J,k,bi,bj)
178     c & + maskC(i,j,Nr,bi,bj)*salt(I,J,Nr,bi,bj)*drF(Nr)*recip_drF(k)
179     salt(I,J,Nr,bi,bj) = 0.
180     ENDDO
181     ENDDO
182    
183     C Note the mapping here is only valid for one tile per proc.
184 adcroft 1.2 DO K = 1, Nr
185     DO J = 1, sNy
186     DO I = 1, sNx
187     I2 = (sNx)*(J-1)+I
188     Katm = _KD2KA( K )
189 cnh 1.4 UG1(I2,Katm,myThid) =
190     & 0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
191     VG1(I2,Katm,myThid) =
192     & 0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))
193     C Physics works with temperature - not potential temp.
194     TG1(I2,Katm,myThid) = theta(I,J,K,bi,bj)
195     & / ((pGround/pSurfs(K))**(RD/CPAIR))
196     c_jmc QG1(I2,Katm,myThid) = salt(I,J,K,bi,bj)
197     QG1(I2,Katm,myThid) = MAX(salt(I,J,K,bi,bj), 0. _d 0)
198     PHIG1(I2,Katm,myThid) = (phiTotal(I,J,K)- ptotalniv5 )
199     & + gravity*Hinitial(k)
200 jmc 1.6 C *NOTE* Fix me for lopped cells <== done !
201     IF (maskC(i,j,k,bi,bj).EQ.1.) THEN
202 cnh 1.4 RHOG1(I2,Katm) = pSurfs(K)/RD/TG1(I2,Katm,myThid)
203 jmc 1.6 ELSE
204 adcroft 1.2 RHOG1(I2,Katm)=0.
205 jmc 1.6 ENDIF
206 adcroft 1.2 ENDDO
207     ENDDO
208     ENDDO
209    
210 cnh 1.4 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
211     c_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
212     DO J = 1, sNy
213     DO I = 1, sNx
214     I2 = I+(J-1)*sNx
215     K = k_surf(i,j,bi,bj)
216     Vsurfsq(I2,myThid) = 0.5 * (
217     & uVel(I,J,K,bi,bj)*uVel(I,J,K,bi,bj)
218     & + uVel(I+1,J,K,bi,bj)*uVel(I+1,J,K,bi,bj)
219     & + vVel(I,J,K,bi,bj)*vVel(I,J,K,bi,bj)
220     & + vVel(I,J+1,K,bi,bj)*vVel(I,J+1,K,bi,bj)
221     & )
222     #ifdef OLD_AIM_GRIG_MAPPING
223     c - to reproduce old results :
224     Katm = _KD2KA( K )
225     Vsurfsq(I2,myThid) =
226     & UG1(I2,Katm,myThid)*UG1(I2,Katm,myThid)
227     & + VG1(I2,Katm,myThid)*VG1(I2,Katm,myThid)
228     #endif /* OLD_AIM_GRIG_MAPPING */
229     ENDDO
230     ENDDO
231     c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
232    
233 adcroft 1.2 C
234     C Set geopotential surfaces
235     C -------------------------
236     DO J=1,sNy
237     DO I=1,sNx
238     I2 = (sNx)*(J-1)+I
239 cnh 1.4 IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
240     PHI0(I2,myThid) = gravity*Hinitialw(Nlevxy(I2,myThid))
241 adcroft 1.2 ELSE
242 cnh 1.4 PHI0(I2,myThid) = 0.
243 adcroft 1.2 ENDIF
244     ENDDO
245     ENDDO
246 cnh 1.4
247 adcroft 1.2 C
248     C Physics package works with log of surface pressure
249     C Get surface pressure from pbot-dpref/dz*Z'
250     DO J=1,sNy
251     DO I=1,sNx
252     I2 = (sNx)*(J-1)+I
253 cnh 1.4 IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
254     PNLEVW(I2,myThid) = PsurfW(Nlevxy(I2,myThid))/pGround
255 adcroft 1.2 ELSE
256     C Dummy value for land
257 cnh 1.4 PNLEVW(I2,myThid) = PsurfW(1)/pGround
258 adcroft 1.2 ENDIF
259 cnh 1.4 PSLG1(I2,myThid) = 0.
260 adcroft 1.2 ENDDO
261     ENDDO
262     cch write(0,*) '(PNLEVW(I2),I2=257,384)'
263     cch write(0,*) (PNLEVW(I2),I2=257,384)
264     C
265     C
266     C Physics package needs to know time of year as a fraction
267     tYear = currentTime/(86400.*360.) -
268     & FLOAT(INT(currentTime/(86400.*360.)))
269 cnh 1.4
270 adcroft 1.2 C
271     C Load external data needed by physics package
272     C 1. Albedo
273     C 2. Soil moisture
274     C 3. Surface temperatures
275     C 4. Snow depth - assume no snow for now
276     C 5. Sea ice - assume no sea ice for now
277     C 6. Land sea mask - infer from exact zeros in soil moisture dataset
278     C 7. Surface geopotential - to be done when orography is in
279     C dynamical kernel. Assume 0. for now.
280     mnthIndex = INT(tYear*12.)+1
281 cnh 1.4 C_cnh01 IF ( mnthIndex .NE. prevMnthIndex .OR.
282     C_cnh01 & FirstCall ) THEN
283     C_cnh01 prevMnthIndex = mnthIndex
284 adcroft 1.2 C Read in surface albedo data (input is in % 0-100 )
285     C scale to give fraction between 0-1 for Francos package.
286 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'salb.',mnthNam(mnthIndex),'.sun.b'
287     C OPEN(1,FILE=fNam(1:14),STATUS='old',FORM='unformatted')
288     C READ(1) tmp4
289     C CLOSE(1)
290 cnh 1.5 C DO J=1,nYio
291     C DO I=1,nXio
292     C tmp4(I,J) = aim_albedo(I,J)/100.
293     C ENDDO
294     C ENDDO
295 adcroft 1.2 DO J=1,sNy
296     DO I=1,sNx
297     I2 = (sNx)*(J-1)+I
298 cnh 1.4 alb0(I2,myThid) = 0.
299 cnh 1.5 alb0(I2,myThid) = aim_albedo(I,J,bi,bj)/100.
300 adcroft 1.2 ENDDO
301     ENDDO
302     C Read in surface temperature data (input is in absolute temperature)
303 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'tsurf.',mnthNam(mnthIndex),'.sun.b'
304     C OPEN(1,FILE=fNam(1:15),STATUS='old',FORM='unformatted')
305     C READ(1) tmp4
306     C CLOSE(1)
307 adcroft 1.2 DO J=1,sNy
308     DO I=1,sNx
309     I2 = (sNx)*(J-1)+I
310 cnh 1.4 sst1(I2,myThid) = 300.
311     stl1(I2,myThid) = 300.
312 cnh 1.5 sst1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
313     stl1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
314 adcroft 1.2 ENDDO
315     ENDDO
316     C
317     C Read in soil moisture data (input is in cm in bucket of depth 20cm. )
318     C??? NOT CLEAR scale for bucket depth of 75mm which is what Franco uses.
319 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'smoist.',mnthNam(mnthIndex),'.sun.b'
320     C OPEN(1,FILE=fNam(1:16),STATUS='old',FORM='unformatted')
321     C READ(1) tmp4
322     C CLOSE(1)
323     C WRITE(0,*) ' Read file ', fNam(1:16), IG0, JG0
324 adcroft 1.2 cdj tmp4 = (tmp4*7.5/20.)*10.
325     DO J=1,sNy
326     DO I=1,sNx
327     I2 = (sNx)*(J-1)+I
328 cnh 1.4 soilq1(I2,myThid) = 0.
329 cnh 1.5 soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)/20.
330 adcroft 1.2 ENDDO
331     ENDDO
332 cnh 1.4 C_cnh01 ENDIF
333 adcroft 1.2 C
334 cnh 1.4 C_cnh01 IF ( FirstCall ) THEN
335 adcroft 1.2 C Set snow depth, sea ice to zero for now
336 cnh 1.4 C Land-sea mask ( figure this out from where
337     C soil moisture is exactly zero ).
338 adcroft 1.2 DO J=1,sNy
339     DO I=1,sNx
340     I2 = (sNx)*(J-1)+I
341 cnh 1.4 fMask1(I2,myThid) = 1.
342     IF ( soilq1(I2,myThid) .EQ. 0. ) fMask1(I2,myThid) = 0.
343     oice1(I2,myThid) = 0.
344     snow1(I2,myThid) = 0.
345 adcroft 1.2 ENDDO
346     ENDDO
347     C open(77,file='lsmask',form='unformatted')
348     C write(77) fmask1
349     C close(77)
350 cnh 1.4 C_cnh01 ENDIF
351 adcroft 1.2 C
352     C Addition may 15 . Reset humidity to 0. if negative
353     C ---------------------------------------------------
354 cnh 1.4 Caja DO K=1,Nr
355     Caja DO J=1-OLy,sNy+OLy
356     Caja DO I=1-Olx,sNx+OLx
357     Caja IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
358     Caja salt(i,j,k,bi,bj) = 0.
359     Caja ENDIF
360     Caja ENDDO
361     Caja ENDDO
362     Caja ENDDO
363 cnh 1.5
364 cnh 1.4
365     CALL PDRIVER( tYear, myThid )
366 adcroft 1.2
367 jmc 1.3 #ifdef ALLOW_TIMEAVE
368 adcroft 1.2 C Calculate diagnostics for AIM
369     CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )
370 jmc 1.3 #endif /* ALLOW_TIMEAVE */
371 adcroft 1.2 C
372     FirstCall = .FALSE.
373 cnh 1.4
374     CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )
375 adcroft 1.2 C
376     #endif /* ALLOW_AIM */
377    
378     RETURN
379     END

  ViewVC Help
Powered by ViewVC 1.1.22