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

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

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

Revision 1.13 - (show annotations) (download)
Mon Sep 29 19:24:31 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint53f_post, checkpoint54a_pre, checkpoint55c_post, checkpoint53b_pre, checkpoint51l_post, checkpoint51j_post, checkpoint57m_post, checkpoint52l_pre, checkpoint52e_pre, hrcube4, hrcube5, checkpoint57g_pre, checkpoint52j_post, checkpoint51o_pre, checkpoint57f_post, checkpoint52e_post, checkpoint51n_pre, checkpoint57j_post, checkpoint57b_post, checkpoint52d_pre, checkpoint53c_post, checkpoint53d_post, checkpoint57f_pre, checkpoint55d_pre, checkpoint57g_post, checkpoint57c_pre, checkpoint55j_post, checkpoint56b_post, checkpoint57h_pre, checkpoint52j_pre, checkpoint54a_post, branch-netcdf, checkpoint55h_post, checkpoint51r_post, checkpoint52b_pre, checkpoint52n_post, checkpoint54b_post, checkpoint51i_post, checkpoint57e_post, checkpoint54d_post, checkpoint56c_post, checkpoint54e_post, checkpoint55b_post, checkpoint51l_pre, checkpoint52m_post, checkpoint55, checkpoint53a_post, checkpoint55a_post, checkpoint57a_post, checkpoint56, checkpoint57o_post, checkpoint55g_post, checkpoint57h_done, checkpoint51o_post, checkpoint57k_post, checkpoint57d_post, checkpoint55f_post, checkpoint57i_post, checkpoint51q_post, checkpoint52l_post, checkpoint52k_post, checkpoint57h_post, checkpoint57a_pre, checkpoint54, checkpoint57, checkpoint53b_post, checkpoint53, checkpoint52, checkpoint52d_post, checkpoint52a_post, checkpoint52b_post, checkpoint53g_post, checkpoint52f_post, checkpoint57n_post, checkpoint52c_post, checkpoint57p_post, checkpoint51h_pre, checkpoint51g_post, ecco_c52_e35, checkpoint54f_post, checkpoint51f_post, eckpoint57e_pre, checkpoint57c_post, checkpoint52a_pre, checkpoint51m_post, checkpoint51t_post, checkpoint53d_pre, checkpoint55e_post, checkpoint54c_post, checkpoint52i_post, checkpoint51p_post, checkpoint51n_post, checkpoint55i_post, checkpoint51i_pre, checkpoint57l_post, checkpoint52i_pre, checkpoint51u_post, checkpoint52h_pre, checkpoint52f_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint56a_post, checkpoint51s_post, checkpoint55d_post
Branch point for: netcdf-sm0, branch-nonh, tg2-branch, checkpoint51n_branch
Changes since 1.12: +2 -2 lines 
 o convert all comments with single quotes (such as "can't", "don't", etc.)
     to unabbreviated form (eg. "do not") since these unmatched quotes
     tend to break the Fortran parser used to generate the HTML-ified
     code browser (see: http://mitgcm.org/dev_docs/code_reference/)
1 C $Header: /u/u3/gcmpack/MITgcm/pkg/aim/aim_do_atmos_physics.F,v 1.12 2002/11/22 03:04:44 jmc Exp $
2 C $Name: $
3
4 #include "AIM_OPTIONS.h"
5
6 SUBROUTINE AIM_DO_PHYSICS( bi, bj, currentTime, myIter, myThid )
7
8 C /==================================================================\
9 C | S/R AIM_DO_ATMOS_PHYSICS |
10 C |==================================================================|
11 C | Interface interface between atmospheric physics package and the |
12 C | dynamical model. |
13 C | Routine calls physics pacakge after mapping model variables to |
14 C | the package grid. Package should derive and set tendency terms |
15 C | which can be included as external forcing terms in the dynamical |
16 C | tendency routines. Packages should communicate this information |
17 C | through common blocks. |
18 C \==================================================================/
19 IMPLICIT NONE
20
21 C -------------- Global variables ------------------------------------
22 C-- size for MITgcm & Physics package :
23 #include "AIM_SIZE.h"
24
25 C-- MITgcm
26 #include "EEPARAMS.h"
27 #include "PARAMS.h"
28 #include "DYNVARS.h"
29 #include "GRID.h"
30 #include "SURFACE.h"
31
32 C-- Physics package
33 #include "AIM_FFIELDS.h"
34 #include "AIM_GRID.h"
35 #include "com_physvar.h"
36 #include "com_forcing1.h"
37
38 C -------------- Routine arguments -----------------------------------
39 c _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
40 _RL currentTime
41 INTEGER myIter, myThid
42 INTEGER bi, bj
43
44 #ifdef ALLOW_AIM
45 C -------------- Local variables -------------------------------------
46 C I,J,K,I2 - Loop counters
47 C tYear - Fraction into year
48 C hInital - Initial height of pressure surfaces (m)
49 C pSurfs - Pressure surfaces (Pa)
50 C Katm - Atmospheric K index
51 INTEGER I
52 INTEGER I2
53 INTEGER J
54 INTEGER K
55 _RL tYear
56 _RL hInitial(Nr)
57 _RL hInitialW(Nr)
58 DATA hInitial / 17338.0,10090.02,5296.88,2038.54,418.038/
59 SAVE hInitial
60 DATA hInitialW / 15090.4, 8050.96, 4087.75, 1657.54, 0. /
61 _RL pSurfs(Nr)
62 DATA pSurfs / 75.D2, 250.D2, 500.D2, 775.D2, 950.D2 /
63 SAVE pSurfs
64 _RL pSurfw(Nr)
65 DATA pSurfw / 150.D2, 350.D2, 650.D2, 900.D2, 1000.D2 /
66 SAVE pSurfw
67 _RL RD
68 _RL CPAIR
69 _RL pGround
70 _RL RhoG1(sNx*sNy,Nr)
71 c _RL phiTotal(sNx,sNy,Nr)
72 c _RL phiTCount
73 c _RL phiTSum
74 c _RL ans
75 c _RL pvoltotNiv5
76 c SAVE pvoltotNiv5
77 c _RL ptotalNiv5
78 INTEGER Katm
79
80 pGround = 1.D5
81 CPAIR = 1004
82 RD = 287
83
84 CALL AIM_DYN2AIM( bi, bj, currentTime, myThid )
85
86 C Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
87 C Internal index mapping is linear in X and Y with a second
88 C dimension for the vertical.
89
90 C Adjustment for heave due to mean heating/cooling
91 C ( I do not think the old formula was strictly "correct" for orography
92 C but I have implemented it as was for now. As implemented
93 C the mean heave of the bottom (K=Nr) level is calculated rather than
94 C the mean heave of the base of the atmosphere. )
95 c phiTCount = 0.
96 c phiTSum = 0.
97 c DO K=1,Nr
98 c DO J=1,sNy
99 c DO I=1,sNx
100 c phiTotal(I,J,K) = etaN(i,j,bi,bj)
101 c phiTCount = phiTCount + hFacC(i,j,Nr,bi,bj)
102 c ENDDO
103 c ENDDO
104 c ENDDO
105 c DO K=1,Nr
106 c DO J=1,sNy
107 c DO I=1,sNx
108 c phiTotal(I,J,K) = phiTotal(I,J,K) +
109 c & recip_rhoConst*(phi_hyd(i,j,k))
110 c ENDDO
111 c ENDDO
112 c ENDDO
113 c DO J=1,sNy
114 c DO I=1,sNx
115 c phiTSum = phiTSum + phiTotal(I,J,Nr)
116 c ENDDO
117 c ENDDO
118 c ans = phiTCount
119 C _GLOBAL_SUM_R8( phiTCount, myThid )
120 c phiTcount = ans
121 c ans = phiTSum
122 C _GLOBAL_SUM_R8( phiTSum, myThid )
123 c phiTSum = ans
124 C ptotalniv5=phiTSum/phiTCount
125 c ptotalniv5=0.
126
127 #ifndef OLD_AIM_INTERFACE
128 C_jmc: Because AIM physics LSC is not applied in the stratosphere (top level),
129 C ==> move water wapor from the stratos to the surface level.
130 DO J = 1-Oly, sNy+Oly
131 DO I = 1-Olx, sNx+Olx
132 k = ksurfC(i,j,bi,bj)
133 IF (k.LE.Nr)
134 & salt(I,J,k,bi,bj) = salt(I,J,k,bi,bj)
135 & + salt(I,J,Nr,bi,bj)*drF(Nr)*recip_drF(k)
136 salt(I,J,Nr,bi,bj) = 0.
137 ENDDO
138 ENDDO
139 #endif /* OLD_AIM_INTERFACE */
140
141 C Note the mapping here is only valid for one tile per proc.
142 DO K = 1, Nr
143 DO J = 1, sNy
144 DO I = 1, sNx
145 I2 = (sNx)*(J-1)+I
146 Katm = _KD2KA( K )
147 C - to reproduce old results (coupled run, summer 2000) :
148 UG1(I2,Katm,myThid) = uVel(I,J,K,bi,bj)
149 VG1(I2,Katm,myThid) = vVel(I,J,K,bi,bj)
150 C Physics works with temperature - not potential temp.
151 TG1(I2,Katm,myThid) = theta(I,J,K,bi,bj)
152 & / ((pGround/pSurfs(Katm))**(RD/CPAIR))
153 #ifdef OLD_AIM_INTERFACE
154 QG1(I2,Katm,myThid) = salt(I,J,K,bi,bj)
155 #else
156 QG1(I2,Katm,myThid) = MAX(salt(I,J,K,bi,bj), 0. _d 0)
157 #endif
158 PHIG1(I2,Katm,myThid) = gravity*Hinitial(Katm)
159 c & +(phiTotal(I,J,K)- ptotalniv5 )
160 IF (maskC(i,j,k,bi,bj).EQ.1.) THEN
161 RHOG1(I2,Katm) = pSurfs(Katm)/RD/TG1(I2,Katm,myThid)
162 ELSE
163 RHOG1(I2,Katm)=0.
164 ENDIF
165 ENDDO
166 ENDDO
167 ENDDO
168
169 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
170 C_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
171 DO J = 1, sNy
172 DO I = 1, sNx
173 I2 = I+(J-1)*sNx
174 #ifdef OLD_AIM_INTERFACE
175 C - to reproduce old results (coupled run, summer 2000) :
176 Vsurfsq(I2,myThid) = 0.
177 IF (NLEVxyU(I2,myThid).GT.0)
178 & Vsurfsq(I2,myThid) = Vsurfsq(I2,myThid)
179 & +UG1(I2,NLEVxyU(I2,myThid),myThid)
180 & *UG1(I2,NLEVxyU(I2,myThid),myThid)
181 IF (NLEVxyV(I2,myThid).GT.0)
182 & Vsurfsq(I2,myThid) = Vsurfsq(I2,myThid)
183 & +VG1(I2,NLEVxyV(I2,myThid),myThid)
184 & *VG1(I2,NLEVxyV(I2,myThid),myThid)
185 #else /* OLD_AIM_INTERFACE */
186 K = ksurfC(i,j,bi,bj)
187 IF (K.LE.Nr) THEN
188 Vsurfsq(I2,myThid) = 0.5 * (
189 & uVel(I,J,K,bi,bj)*uVel(I,J,K,bi,bj)
190 & + uVel(I+1,J,K,bi,bj)*uVel(I+1,J,K,bi,bj)
191 & + vVel(I,J,K,bi,bj)*vVel(I,J,K,bi,bj)
192 & + vVel(I,J+1,K,bi,bj)*vVel(I,J+1,K,bi,bj)
193 & )
194 ELSE
195 Vsurfsq(I2,myThid) = 0.
196 ENDIF
197 #endif /* OLD_AIM_INTERFACE */
198 ENDDO
199 ENDDO
200 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
201
202 C
203 C Set geopotential surfaces
204 C -------------------------
205 DO J=1,sNy
206 DO I=1,sNx
207 I2 = (sNx)*(J-1)+I
208 IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
209 PHI0(I2,myThid) = gravity*Hinitialw(Nlevxy(I2,myThid))
210 ELSE
211 PHI0(I2,myThid) = 0.
212 ENDIF
213 ENDDO
214 ENDDO
215
216 C
217 C Physics package works with log of surface pressure
218 C Get surface pressure from pbot-dpref/dz*Z'
219 DO J=1,sNy
220 DO I=1,sNx
221 I2 = (sNx)*(J-1)+I
222 IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
223 PNLEVW(I2,myThid) = PsurfW(Nlevxy(I2,myThid))/pGround
224 ELSE
225 C Dummy value for land
226 PNLEVW(I2,myThid) = PsurfW(Nr)/pGround
227 ENDIF
228 PSLG1(I2,myThid) = 0.
229 ENDDO
230 ENDDO
231 C
232 C
233 C Physics package needs to know time of year as a fraction
234 tYear = currentTime/(86400.*360.) -
235 & FLOAT(INT(currentTime/(86400.*360.)))
236
237 C
238 C Load external data needed by physics package
239 C 1. Albedo (between 0-1)
240 C 2. Soil moisture (between 0-1)
241 C 3. Surface temperatures (in situ Temp. [K])
242 C 4. Snow depth - assume no snow for now
243 C 5. Sea ice - assume no sea ice for now
244 C 6. Land sea mask - infer from exact zeros in soil moisture dataset
245 C 7. Surface geopotential - to be done when orography is in
246 C dynamical kernel. Assume 0. for now.
247
248 C Load in surface albedo data (in [0,1]) from aim_albedo to alb0 :
249 DO J=1,sNy
250 DO I=1,sNx
251 I2 = (sNx)*(J-1)+I
252 alb0(I2,myThid) = 0.
253 alb0(I2,myThid) = aim_albedo(I,J,bi,bj)
254 ENDDO
255 ENDDO
256 C Load in surface temperature data from aim_surfTemp to stl1 & sst1 :
257 DO J=1,sNy
258 DO I=1,sNx
259 I2 = (sNx)*(J-1)+I
260 sst1(I2,myThid) = 300.
261 stl1(I2,myThid) = 300.
262 sst1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
263 stl1(I2,myThid) = aim_surfTemp(I,J,bi,bj)
264 ENDDO
265 ENDDO
266 C
267 C Load in soil moisture data (in [0,1]) from aim_soilMoisture to soilq1 :
268 C??? NOT CLEAR scale for bucket depth of 75mm which is what Franco uses.
269 DO J=1,sNy
270 DO I=1,sNx
271 I2 = (sNx)*(J-1)+I
272 soilq1(I2,myThid) = 0.
273 soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)
274 ENDDO
275 ENDDO
276 C Set snow depth, sea ice to zero for now
277 C Land-sea mask ( figure this out from where
278 C soil moisture is exactly zero ).
279 DO J=1,sNy
280 DO I=1,sNx
281 I2 = (sNx)*(J-1)+I
282 fMask1(I2,myThid) = 1.
283 IF ( soilq1(I2,myThid) .EQ. 0. ) fMask1(I2,myThid) = 0.
284 oice1(I2,myThid) = 0.
285 snow1(I2,myThid) = 0.
286 ENDDO
287 ENDDO
288 C open(77,file='lsmask',form='unformatted')
289 C write(77) fmask1
290 C close(77)
291
292 C Addition may 15 . Reset humidity to 0. if negative
293 C ---------------------------------------------------
294 #ifdef OLD_AIM_INTERFACE
295 DO K=1,Nr
296 DO J=1-OLy,sNy+OLy
297 DO I=1-Olx,sNx+OLx
298 IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
299 salt(i,j,k,bi,bj) = 0.
300 ENDIF
301 ENDDO
302 ENDDO
303 ENDDO
304 #endif /* OLD_AIM_INTERFACE */
305
306 CALL PDRIVER( tYear, myThid )
307
308 #ifdef ALLOW_TIMEAVE
309 C Calculate diagnostics for AIM
310 CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )
311 #endif /* ALLOW_TIMEAVE */
312
313 CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )
314
315 #endif /* ALLOW_AIM */
316
317 RETURN
318 END
  ViewVC Help
Powered by ViewVC 1.1.22