/[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.2.1 - (hide annotations) (download)
Tue Feb 26 16:04:48 2002 UTC (22 years, 2 months ago) by adcroft
Branch: release1
CVS Tags: release1_p12, release1_p10, release1_p11, release1_p16, release1_p15, release1_p13_pre, release1_p14, release1_p12_pre, release1_p13, release1_p17, release1_p9, release1_p8, release1_p2, release1_p3, release1_p4, release1_p6, release1_p1, release1_p5, release1_p7, release1_chkpt44d_post
Branch point for: release1_50yr
Changes since 1.8: +12 -9 lines
Merging changes on MAIN between checkpoint43 and checkpoint43a-release1mods
Command: cvs -q update -jcheckpoint43 -jcheckpoint43a-release1mods -d -P

These changes are most of the changes between c43 and c44 except those
that occured after "12:45 11 Jan 2002". As far as I can tell it is
checkpoint43 with the following mods:

  o fix bug in mom_vi_del2uv
  o select when filters are applied ; add options to zonal_filter (data.zonfilt)  o gmredi: fix Pb in the adiabatic form ; add options (.e.g. Bolus advection)
  o update AIM experiments (NCEP input files)
  o improve and extend diagnostics (Monitor, TimeAve with NonLin-FrSurf)
  o added some stuff for AD
  o Jamar wet-points

This update does not contain the following mods that are in checkpoint44

  o bug fix in pkg/generic_advdiff/
    - thread related bug, bi,bj arguments in vertical advection routines
  o some changes to pkg/autodiff, pkg/cost, pkg/exf, pkg/ecco,
    verification/carbon and model/src/ related to adjoint
  o some new Matlab scripts for diagnosing model density
    - utils/matlab/dens_poly3.m and ini_poly3.m

The list of exclusions is accurate based on a "cvs diff". The list of
inclusions is based on the record in doc/tag-index which may not be complete.

1 adcroft 1.8.2.1 C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/aim/Attic/aim_do_atmos_physics.F,v 1.8.2.1 2002/02/26 16:04:48 adcroft 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 adcroft 1.8.2.1 k = ksurfC(i,j,bi,bj)
177     IF (k.LE.Nr)
178     & salt(I,J,k,bi,bj) = salt(I,J,k,bi,bj)
179     & + salt(I,J,Nr,bi,bj)*drF(Nr)*recip_drF(k)
180 cnh 1.4 salt(I,J,Nr,bi,bj) = 0.
181     ENDDO
182     ENDDO
183 jmc 1.8 #endif /* OLD_AIM_INTERFACE */
184 cnh 1.4
185     C Note the mapping here is only valid for one tile per proc.
186 adcroft 1.2 DO K = 1, Nr
187     DO J = 1, sNy
188     DO I = 1, sNx
189     I2 = (sNx)*(J-1)+I
190     Katm = _KD2KA( K )
191 jmc 1.8 C - to reproduce old results (coupled run, summer 2000) :
192     UG1(I2,Katm,myThid) = uVel(I,J,K,bi,bj)
193     VG1(I2,Katm,myThid) = vVel(I,J,K,bi,bj)
194 cnh 1.4 C Physics works with temperature - not potential temp.
195     TG1(I2,Katm,myThid) = theta(I,J,K,bi,bj)
196 jmc 1.8 & / ((pGround/pSurfs(Katm))**(RD/CPAIR))
197     #ifdef OLD_AIM_INTERFACE
198     QG1(I2,Katm,myThid) = salt(I,J,K,bi,bj)
199     #else
200 cnh 1.4 QG1(I2,Katm,myThid) = MAX(salt(I,J,K,bi,bj), 0. _d 0)
201 jmc 1.8 #endif
202 cnh 1.4 PHIG1(I2,Katm,myThid) = (phiTotal(I,J,K)- ptotalniv5 )
203 jmc 1.8 & + gravity*Hinitial(Katm)
204 jmc 1.6 C *NOTE* Fix me for lopped cells <== done !
205     IF (maskC(i,j,k,bi,bj).EQ.1.) THEN
206 jmc 1.8 RHOG1(I2,Katm) = pSurfs(Katm)/RD/TG1(I2,Katm,myThid)
207 jmc 1.6 ELSE
208 adcroft 1.2 RHOG1(I2,Katm)=0.
209 jmc 1.6 ENDIF
210 adcroft 1.2 ENDDO
211     ENDDO
212     ENDDO
213    
214 cnh 1.4 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
215     c_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
216     DO J = 1, sNy
217     DO I = 1, sNx
218 jmc 1.7 I2 = I+(J-1)*sNx
219 jmc 1.8 #ifdef OLD_AIM_INTERFACE
220     C - to reproduce old results (coupled run, summer 2000) :
221     Vsurfsq(I2,myThid) = 0.
222     IF (NLEVxyU(I2,myThid).GT.0)
223     & Vsurfsq(I2,myThid) = Vsurfsq(I2,myThid)
224     & +UG1(I2,NLEVxyU(I2,myThid),myThid)
225     & *UG1(I2,NLEVxyU(I2,myThid),myThid)
226     IF (NLEVxyV(I2,myThid).GT.0)
227     & Vsurfsq(I2,myThid) = Vsurfsq(I2,myThid)
228     & +VG1(I2,NLEVxyV(I2,myThid),myThid)
229     & *VG1(I2,NLEVxyV(I2,myThid),myThid)
230     #else /* OLD_AIM_INTERFACE */
231 jmc 1.7 K = ksurfC(i,j,bi,bj)
232     IF (K.LE.Nr) THEN
233 cnh 1.4 Vsurfsq(I2,myThid) = 0.5 * (
234     & uVel(I,J,K,bi,bj)*uVel(I,J,K,bi,bj)
235     & + uVel(I+1,J,K,bi,bj)*uVel(I+1,J,K,bi,bj)
236     & + vVel(I,J,K,bi,bj)*vVel(I,J,K,bi,bj)
237     & + vVel(I,J+1,K,bi,bj)*vVel(I,J+1,K,bi,bj)
238     & )
239 jmc 1.7 ELSE
240     Vsurfsq(I2,myThid) = 0.
241     ENDIF
242 jmc 1.8 #endif /* OLD_AIM_INTERFACE */
243 cnh 1.4 ENDDO
244     ENDDO
245     c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
246    
247 adcroft 1.2 C
248     C Set geopotential surfaces
249     C -------------------------
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     PHI0(I2,myThid) = gravity*Hinitialw(Nlevxy(I2,myThid))
255 adcroft 1.2 ELSE
256 cnh 1.4 PHI0(I2,myThid) = 0.
257 adcroft 1.2 ENDIF
258     ENDDO
259     ENDDO
260 cnh 1.4
261 adcroft 1.2 C
262     C Physics package works with log of surface pressure
263     C Get surface pressure from pbot-dpref/dz*Z'
264     DO J=1,sNy
265     DO I=1,sNx
266     I2 = (sNx)*(J-1)+I
267 cnh 1.4 IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
268     PNLEVW(I2,myThid) = PsurfW(Nlevxy(I2,myThid))/pGround
269 adcroft 1.2 ELSE
270     C Dummy value for land
271 jmc 1.8 PNLEVW(I2,myThid) = PsurfW(Nr)/pGround
272 adcroft 1.2 ENDIF
273 cnh 1.4 PSLG1(I2,myThid) = 0.
274 adcroft 1.2 ENDDO
275     ENDDO
276     cch write(0,*) '(PNLEVW(I2),I2=257,384)'
277     cch write(0,*) (PNLEVW(I2),I2=257,384)
278     C
279     C
280     C Physics package needs to know time of year as a fraction
281     tYear = currentTime/(86400.*360.) -
282     & FLOAT(INT(currentTime/(86400.*360.)))
283 cnh 1.4
284 adcroft 1.2 C
285     C Load external data needed by physics package
286 adcroft 1.8.2.1 C 1. Albedo (between 0-1)
287     C 2. Soil moisture (between 0-1)
288     C 3. Surface temperatures (in situ Temp. [K])
289 adcroft 1.2 C 4. Snow depth - assume no snow for now
290     C 5. Sea ice - assume no sea ice for now
291     C 6. Land sea mask - infer from exact zeros in soil moisture dataset
292     C 7. Surface geopotential - to be done when orography is in
293     C dynamical kernel. Assume 0. for now.
294     mnthIndex = INT(tYear*12.)+1
295 cnh 1.4 C_cnh01 IF ( mnthIndex .NE. prevMnthIndex .OR.
296     C_cnh01 & FirstCall ) THEN
297     C_cnh01 prevMnthIndex = mnthIndex
298 adcroft 1.2 C Read in surface albedo data (input is in % 0-100 )
299     C scale to give fraction between 0-1 for Francos package.
300 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'salb.',mnthNam(mnthIndex),'.sun.b'
301     C OPEN(1,FILE=fNam(1:14),STATUS='old',FORM='unformatted')
302     C READ(1) tmp4
303     C CLOSE(1)
304 cnh 1.5 C DO J=1,nYio
305     C DO I=1,nXio
306     C tmp4(I,J) = aim_albedo(I,J)/100.
307     C ENDDO
308     C ENDDO
309 adcroft 1.2 DO J=1,sNy
310     DO I=1,sNx
311     I2 = (sNx)*(J-1)+I
312 cnh 1.4 alb0(I2,myThid) = 0.
313 adcroft 1.8.2.1 c alb0(I2,myThid) = aim_albedo(I,J,bi,bj)/100.
314     alb0(I2,myThid) = aim_albedo(I,J,bi,bj)
315 adcroft 1.2 ENDDO
316     ENDDO
317     C Read in surface temperature data (input is in absolute temperature)
318 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'tsurf.',mnthNam(mnthIndex),'.sun.b'
319     C OPEN(1,FILE=fNam(1:15),STATUS='old',FORM='unformatted')
320     C READ(1) tmp4
321     C CLOSE(1)
322 adcroft 1.2 DO J=1,sNy
323     DO I=1,sNx
324     I2 = (sNx)*(J-1)+I
325 cnh 1.4 sst1(I2,myThid) = 300.
326     stl1(I2,myThid) = 300.
327 cnh 1.5 sst1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
328     stl1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
329 adcroft 1.2 ENDDO
330     ENDDO
331     C
332     C Read in soil moisture data (input is in cm in bucket of depth 20cm. )
333     C??? NOT CLEAR scale for bucket depth of 75mm which is what Franco uses.
334 cnh 1.4 C WRITE(fNam,'(A,A,A)' ) 'smoist.',mnthNam(mnthIndex),'.sun.b'
335     C OPEN(1,FILE=fNam(1:16),STATUS='old',FORM='unformatted')
336     C READ(1) tmp4
337     C CLOSE(1)
338     C WRITE(0,*) ' Read file ', fNam(1:16), IG0, JG0
339 adcroft 1.2 cdj tmp4 = (tmp4*7.5/20.)*10.
340     DO J=1,sNy
341     DO I=1,sNx
342     I2 = (sNx)*(J-1)+I
343 cnh 1.4 soilq1(I2,myThid) = 0.
344 adcroft 1.8.2.1 c soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)/20.
345     soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)
346 adcroft 1.2 ENDDO
347     ENDDO
348 cnh 1.4 C_cnh01 ENDIF
349 adcroft 1.2 C
350 cnh 1.4 C_cnh01 IF ( FirstCall ) THEN
351 adcroft 1.2 C Set snow depth, sea ice to zero for now
352 cnh 1.4 C Land-sea mask ( figure this out from where
353     C soil moisture is exactly zero ).
354 adcroft 1.2 DO J=1,sNy
355     DO I=1,sNx
356     I2 = (sNx)*(J-1)+I
357 cnh 1.4 fMask1(I2,myThid) = 1.
358     IF ( soilq1(I2,myThid) .EQ. 0. ) fMask1(I2,myThid) = 0.
359     oice1(I2,myThid) = 0.
360     snow1(I2,myThid) = 0.
361 adcroft 1.2 ENDDO
362     ENDDO
363     C open(77,file='lsmask',form='unformatted')
364     C write(77) fmask1
365     C close(77)
366 cnh 1.4 C_cnh01 ENDIF
367 adcroft 1.2 C
368     C Addition may 15 . Reset humidity to 0. if negative
369     C ---------------------------------------------------
370 jmc 1.8 #ifdef OLD_AIM_INTERFACE
371     DO K=1,Nr
372     DO J=1-OLy,sNy+OLy
373     DO I=1-Olx,sNx+OLx
374     IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
375     salt(i,j,k,bi,bj) = 0.
376     ENDIF
377     ENDDO
378     ENDDO
379     ENDDO
380     #endif /* OLD_AIM_INTERFACE */
381 cnh 1.4
382     CALL PDRIVER( tYear, myThid )
383 adcroft 1.2
384 jmc 1.3 #ifdef ALLOW_TIMEAVE
385 adcroft 1.2 C Calculate diagnostics for AIM
386     CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )
387 jmc 1.3 #endif /* ALLOW_TIMEAVE */
388 adcroft 1.2 C
389     FirstCall = .FALSE.
390 cnh 1.4
391     CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )
392 adcroft 1.2 C
393     #endif /* ALLOW_AIM */
394    
395     RETURN
396     END

  ViewVC Help
Powered by ViewVC 1.1.22