/[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.5 - (hide annotations) (download)
Mon Jun 18 17:39:58 2001 UTC (22 years, 11 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint40pre1, checkpoint40pre2, checkpoint40pre5, checkpoint40pre6, checkpoint40pre4, checkpoint40pre3, checkpoint40pre7
Changes since 1.4: +10 -15 lines
Add to main branch of
  o CS atmos with AIM physics
  o Multi-threaded AIM physics for LatLon and CS tests
  o Tidied up monitor() output

1 cnh 1.4 C $Header: /u/gcmpack/models/MITgcmUV/verification/aim.5l_LatLon/code/Attic/aim_do_atmos_physics.F,v 1.1.2.2 2001/04/17 01:29:15 jmc 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     C *NOTE* Fix me for lopped cells
201 adcroft 1.2 if(hFacC(i,j,k,bi,bj).eq.1.) then
202 cnh 1.4 RHOG1(I2,Katm) = pSurfs(K)/RD/TG1(I2,Katm,myThid)
203 adcroft 1.2 else
204     RHOG1(I2,Katm)=0.
205     endif
206 cnh 1.4 C *NOTE* Fix me for lopped cells
207 adcroft 1.2 ENDDO
208     ENDDO
209     ENDDO
210    
211 cnh 1.4 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
212     c_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
213     DO J = 1, sNy
214     DO I = 1, sNx
215     I2 = I+(J-1)*sNx
216     K = k_surf(i,j,bi,bj)
217     Vsurfsq(I2,myThid) = 0.5 * (
218     & uVel(I,J,K,bi,bj)*uVel(I,J,K,bi,bj)
219     & + uVel(I+1,J,K,bi,bj)*uVel(I+1,J,K,bi,bj)
220     & + vVel(I,J,K,bi,bj)*vVel(I,J,K,bi,bj)
221     & + vVel(I,J+1,K,bi,bj)*vVel(I,J+1,K,bi,bj)
222     & )
223     #ifdef OLD_AIM_GRIG_MAPPING
224     c - to reproduce old results :
225     Katm = _KD2KA( K )
226     Vsurfsq(I2,myThid) =
227     & UG1(I2,Katm,myThid)*UG1(I2,Katm,myThid)
228     & + VG1(I2,Katm,myThid)*VG1(I2,Katm,myThid)
229     #endif /* OLD_AIM_GRIG_MAPPING */
230     ENDDO
231     ENDDO
232     c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
233    
234 adcroft 1.2 C
235     C Set geopotential surfaces
236     C -------------------------
237     DO J=1,sNy
238     DO I=1,sNx
239     I2 = (sNx)*(J-1)+I
240 cnh 1.4 IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
241     PHI0(I2,myThid) = gravity*Hinitialw(Nlevxy(I2,myThid))
242 adcroft 1.2 ELSE
243 cnh 1.4 PHI0(I2,myThid) = 0.
244 adcroft 1.2 ENDIF
245     ENDDO
246     ENDDO
247 cnh 1.4
248 adcroft 1.2 C
249     C Physics package works with log of surface pressure
250     C Get surface pressure from pbot-dpref/dz*Z'
251     DO J=1,sNy
252     DO I=1,sNx
253     I2 = (sNx)*(J-1)+I
254 cnh 1.4 IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
255     PNLEVW(I2,myThid) = PsurfW(Nlevxy(I2,myThid))/pGround
256 adcroft 1.2 ELSE
257     C Dummy value for land
258 cnh 1.4 PNLEVW(I2,myThid) = PsurfW(1)/pGround
259 adcroft 1.2 ENDIF
260 cnh 1.4 PSLG1(I2,myThid) = 0.
261 adcroft 1.2 ENDDO
262     ENDDO
263     cch write(0,*) '(PNLEVW(I2),I2=257,384)'
264     cch write(0,*) (PNLEVW(I2),I2=257,384)
265     C
266     C
267     C Physics package needs to know time of year as a fraction
268     tYear = currentTime/(86400.*360.) -
269     & FLOAT(INT(currentTime/(86400.*360.)))
270 cnh 1.4
271 adcroft 1.2 C
272     C Load external data needed by physics package
273     C 1. Albedo
274     C 2. Soil moisture
275     C 3. Surface temperatures
276     C 4. Snow depth - assume no snow for now
277     C 5. Sea ice - assume no sea ice for now
278     C 6. Land sea mask - infer from exact zeros in soil moisture dataset
279     C 7. Surface geopotential - to be done when orography is in
280     C dynamical kernel. Assume 0. for now.
281     mnthIndex = INT(tYear*12.)+1
282 cnh 1.4 C_cnh01 IF ( mnthIndex .NE. prevMnthIndex .OR.
283     C_cnh01 & FirstCall ) THEN
284     C_cnh01 prevMnthIndex = mnthIndex
285 adcroft 1.2 C Read in surface albedo data (input is in % 0-100 )
286     C scale to give fraction between 0-1 for Francos package.
287 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'salb.',mnthNam(mnthIndex),'.sun.b'
288     C OPEN(1,FILE=fNam(1:14),STATUS='old',FORM='unformatted')
289     C READ(1) tmp4
290     C CLOSE(1)
291 cnh 1.5 C DO J=1,nYio
292     C DO I=1,nXio
293     C tmp4(I,J) = aim_albedo(I,J)/100.
294     C ENDDO
295     C ENDDO
296 adcroft 1.2 DO J=1,sNy
297     DO I=1,sNx
298     I2 = (sNx)*(J-1)+I
299 cnh 1.4 alb0(I2,myThid) = 0.
300 cnh 1.5 alb0(I2,myThid) = aim_albedo(I,J,bi,bj)/100.
301 adcroft 1.2 ENDDO
302     ENDDO
303     C Read in surface temperature data (input is in absolute temperature)
304 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'tsurf.',mnthNam(mnthIndex),'.sun.b'
305     C OPEN(1,FILE=fNam(1:15),STATUS='old',FORM='unformatted')
306     C READ(1) tmp4
307     C CLOSE(1)
308 adcroft 1.2 DO J=1,sNy
309     DO I=1,sNx
310     I2 = (sNx)*(J-1)+I
311 cnh 1.4 sst1(I2,myThid) = 300.
312     stl1(I2,myThid) = 300.
313 cnh 1.5 sst1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
314     stl1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
315 adcroft 1.2 ENDDO
316     ENDDO
317     C
318     C Read in soil moisture data (input is in cm in bucket of depth 20cm. )
319     C??? NOT CLEAR scale for bucket depth of 75mm which is what Franco uses.
320 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'smoist.',mnthNam(mnthIndex),'.sun.b'
321     C OPEN(1,FILE=fNam(1:16),STATUS='old',FORM='unformatted')
322     C READ(1) tmp4
323     C CLOSE(1)
324     C WRITE(0,*) ' Read file ', fNam(1:16), IG0, JG0
325 adcroft 1.2 cdj tmp4 = (tmp4*7.5/20.)*10.
326     DO J=1,sNy
327     DO I=1,sNx
328     I2 = (sNx)*(J-1)+I
329 cnh 1.4 soilq1(I2,myThid) = 0.
330 cnh 1.5 soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)/20.
331 adcroft 1.2 ENDDO
332     ENDDO
333 cnh 1.4 C_cnh01 ENDIF
334 adcroft 1.2 C
335 cnh 1.4 C_cnh01 IF ( FirstCall ) THEN
336 adcroft 1.2 C Set snow depth, sea ice to zero for now
337 cnh 1.4 C Land-sea mask ( figure this out from where
338     C soil moisture is exactly zero ).
339 adcroft 1.2 DO J=1,sNy
340     DO I=1,sNx
341     I2 = (sNx)*(J-1)+I
342 cnh 1.4 fMask1(I2,myThid) = 1.
343     IF ( soilq1(I2,myThid) .EQ. 0. ) fMask1(I2,myThid) = 0.
344     oice1(I2,myThid) = 0.
345     snow1(I2,myThid) = 0.
346 adcroft 1.2 ENDDO
347     ENDDO
348     C open(77,file='lsmask',form='unformatted')
349     C write(77) fmask1
350     C close(77)
351 cnh 1.4 C_cnh01 ENDIF
352 adcroft 1.2 C
353     C Addition may 15 . Reset humidity to 0. if negative
354     C ---------------------------------------------------
355 cnh 1.4 Caja DO K=1,Nr
356     Caja DO J=1-OLy,sNy+OLy
357     Caja DO I=1-Olx,sNx+OLx
358     Caja IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
359     Caja salt(i,j,k,bi,bj) = 0.
360     Caja ENDIF
361     Caja ENDDO
362     Caja ENDDO
363     Caja ENDDO
364 cnh 1.5
365 cnh 1.4
366     CALL PDRIVER( tYear, myThid )
367 adcroft 1.2
368 jmc 1.3 #ifdef ALLOW_TIMEAVE
369 adcroft 1.2 C Calculate diagnostics for AIM
370     CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )
371 jmc 1.3 #endif /* ALLOW_TIMEAVE */
372 adcroft 1.2 C
373     FirstCall = .FALSE.
374 cnh 1.4
375     CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )
376 adcroft 1.2 C
377     #endif /* ALLOW_AIM */
378    
379     RETURN
380     END

  ViewVC Help
Powered by ViewVC 1.1.22