/[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.8 - (hide annotations) (download)
Tue Sep 25 19:55:34 2001 UTC (22 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: release1_b1, checkpoint42, checkpoint43, ecco-branch-mod1, release1_beta1
Branch point for: ecco-branch, release1_coupled, release1
Changes since 1.7: +42 -32 lines
add a CPP option to turn back to old AIM Interface.
diagnostic of surface stress consistent with dynamical effect.
bug with Katm fixed.

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

  ViewVC Help
Powered by ViewVC 1.1.22