/[MITgcm]/MITgcm/pkg/aim_v23/phy_driver.F
ViewVC logotype

Contents of /MITgcm/pkg/aim_v23/phy_driver.F

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


Revision 1.10 - (show annotations) (download)
Thu Jan 11 01:58:50 2018 UTC (6 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, HEAD
Changes since 1.9: +15 -5 lines
set LW absorption in CO2 band (prescribed or computed from pCO2)
in S/R PHY_DRIVER and pass it as argument to S/R RADSW

1 C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/phy_driver.F,v 1.9 2010/10/26 20:59:53 dfer Exp $
2 C $Name: $
3
4 #include "AIM_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: PHY_DRIVER
8 C !INTERFACE:
9 SUBROUTINE PHY_DRIVER( tYear, usePkgDiag,
10 I bi, bj, myTime, myIter, myThid )
11
12 C !DESCRIPTION: \bv
13 C------------------------
14 C-- SUBROUTINE PHYDRIVER (tYear, myTime, bi, bj, myThid )
15 C-- Purpose: stand-alone driver for physical parametrization routines
16 C-- Input : TYEAR : fraction of year (0 = 1jan.00, 1 = 31dec.24)
17 C-- grid-point model fields in common block: PHYGR1
18 C-- forcing fields in common blocks : LSMASK, FORFIX, FORCIN
19 C-- Output : Diagnosed upper-air variables in common block: PHYGR2
20 C-- Diagnosed surface variables in common block: PHYGR3
21 C-- Physical param. tendencies in common block: PHYTEN
22 C-- Surface and upper boundary fluxes in common block: FLUXES
23 C-------
24 C Note: tendencies are not /dpFac here but later in AIM_AIM2DYN
25 C-------
26 C from SPEDDY code: (part of original code left with c_FM)
27 C * S/R PHYPAR : except interp. dynamical Var. from Spectral of grid point
28 C here dynamical var. are loaded within S/R AIM_DYN2AIM.
29 C * S/R FORDATE: only the CALL SOL_OZ (done once / day in SPEEDY)
30 C------------------------
31 C \ev
32
33 C !USES:
34 IMPLICIT NONE
35
36 C == Global variables ===
37
38 C-- size for MITgcm & Physics package :
39 #include "AIM_SIZE.h"
40 #include "EEPARAMS.h"
41
42 C-- Physics package
43 #include "AIM_PARAMS.h"
44 #include "AIM_GRID.h"
45 #include "AIM_CO2.h"
46
47 C Constants + functions of sigma and latitude
48 #include "com_physcon.h"
49
50 C Model variables, tendencies and fluxes on gaussian grid
51 #include "com_physvar.h"
52
53 C Surface forcing fields (time-inv. or functions of seasonal cycle)
54 #include "com_forcing.h"
55
56 C Constants for forcing fields:
57 #include "com_forcon.h"
58
59 C Radiation scheme variables
60 #include "com_radvar.h"
61
62 C Radiation constants
63 #include "com_radcon.h"
64
65 C Logical flags
66 c_FM include "com_lflags.h"
67
68 C !INPUT/OUTPUT PARAMETERS:
69 C == Routine arguments ==
70 C tYear :: Fraction into year
71 C usePkgDiag :: logical flag, true if using Diagnostics PKG
72 C bi, bj :: Tile index
73 C myTime :: Current time of simulation ( s )
74 C myIter :: Current iteration number in simulation
75 C myThid :: Number of this instance of the routine
76 _RL tYear
77 LOGICAL usePkgDiag
78 INTEGER bi,bj
79 _RL myTime
80 INTEGER myIter, myThid
81 CEOP
82
83 #ifdef ALLOW_AIM
84 C !FUNCTIONS:
85 C !LOCAL VARIABLES:
86 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
87 C-- Local Variables originally (Speedy) in common bloc (com_physvar.h):
88 C TG1 :: absolute temperature
89 C QG1 :: specific humidity (g/kg)
90 C VsurfSq :: Square of surface wind speed (grid position = as T,Q)
91 C SE :: dry static energy <- replaced by Pot.Temp.
92 C QSAT :: saturation specific humidity (g/kg)
93 C PSG :: surface pressure (normalized)
94 _RL TG1 (NGP,NLEV)
95 _RL QG1 (NGP,NLEV)
96 _RL VsurfSq(NGP)
97 _RL SE (NGP,NLEV)
98 _RL QSAT (NGP,NLEV)
99 _RL PSG (NGP)
100 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
101 C-- Local variables:
102 C absLW_CO2 :: LW absorbtion in CO2 band (uniform value)
103 C kGround :: Ground level index (2-dim)
104 C dpFac :: cell delta_P fraction (3-dim)
105 C dTskin :: temp. correction for daily-cycle heating [K]
106 C T1s :: near-surface air temperature (from Pot.Temp)
107 C DENVV :: surface flux (sens,lat.) coeff. (=Rho*|V|) [kg/m2/s]
108 C Shf0 :: sensible heat flux over freezing surf.
109 C dShf :: sensible heat flux derivative relative to surf. temp
110 C Evp0 :: evaporation computed over freezing surface (Ts=0.oC)
111 C dEvp :: evaporation derivative relative to surf. temp
112 C Slr0 :: upward long wave radiation over freezing surf.
113 C dSlr :: upward long wave rad. derivative relative to surf. temp
114 C sFlx :: net surface flux (+=down) function of surf. temp Ts:
115 C 0: Flux(Ts=0.oC) ; 1: Flux(Ts^n) ; 2: d.Flux/d.Ts(Ts^n)
116 LOGICAL LRADSW
117 INTEGER ICLTOP(NGP)
118 INTEGER kGround(NGP)
119 _RL absLW_CO2
120 _RL dpFac(NGP,NLEV)
121 c_FM REAL RPS(NGP), ST4S(NGP)
122 _RL ST4S(NGP)
123 _RL PSG_1(NGP), RPS_1
124 _RL dTskin(NGP), T1s(NGP), DENVV(NGP)
125 _RL Shf0(NGP), dShf(NGP), Evp0(NGP), dEvp(NGP)
126 _RL Slr0(NGP), dSlr(NGP), sFlx(NGP,0:2)
127 _RL UPSWG(NGP)
128
129 INTEGER J, K
130
131 #ifdef ALLOW_CLR_SKY_DIAG
132 _RL dummyR(NGP)
133 INTEGER dummyI(NGP)
134 #endif
135 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
136
137 C-- 1. Compute grid-point fields
138
139 C- 1.1 Convert model spectral variables to grid-point variables
140
141 CALL AIM_DYN2AIM(
142 O TG1, QG1, SE, VsurfSq, PSG, dpFac, kGround,
143 I bi, bj, myTime, myIter, myThid )
144
145 C- 1.2 Compute thermodynamic variables
146
147 C- 1.2.a Surface pressure (ps), 1/ps and surface temperature
148 RPS_1 = 1. _d 0
149 DO J=1,NGP
150 PSG_1(J)=1. _d 0
151 c_FM PSG(J)=EXP(PSLG1(J))
152 c_FM RPS(J)=1./PSG(J)
153 ENDDO
154
155 C 1.2.b Dry static energy
156 C <= replaced by Pot.Temp in aim_dyn2aim
157 c DO K=1,NLEV
158 c DO J=1,NGP
159 c_FM SE(J,K)=CP*TG1(J,K)+PHIG1(J,K)
160 c ENDDO
161 c ENDDO
162
163 C 1.2.c Relative humidity and saturation spec. humidity
164
165 DO K=1,NLEV
166 c_FM CALL SHTORH (1,NGP,TG1(1,K),PSG,SIG(K),QG1(1,K),
167 c_FM & RH(1,K),QSAT(1,K))
168 CALL SHTORH (1,NGP,TG1(1,K),PSG_1,SIG(K),QG1(1,K),
169 O RH(1,K,myThid),QSAT(1,K),
170 I myThid)
171 ENDDO
172
173 C-- 2. Precipitation
174
175 C 2.1 Deep convection
176
177 c_FM CALL CONVMF (PSG,SE,QG1,QSAT,
178 c_FM & ICLTOP,CBMF,PRECNV,TT_CNV,QT_CNV)
179 CALL CONVMF (PSG,dpFac,SE,QG1,QSAT,
180 O ICLTOP,CBMF(1,myThid),PRECNV(1,myThid),
181 O TT_CNV(1,1,myThid),QT_CNV(1,1,myThid),
182 I kGround,bi,bj,myThid)
183
184 DO K=2,NLEV
185 DO J=1,NGP
186 TT_CNV(J,K,myThid)=TT_CNV(J,K,myThid)*RPS_1*GRDSCP(K)
187 QT_CNV(J,K,myThid)=QT_CNV(J,K,myThid)*RPS_1*GRDSIG(K)
188 ENDDO
189 ENDDO
190
191 C 2.2 Large-scale condensation
192
193 c_FM CALL LSCOND (PSG,QG1,QSAT,
194 c_FM & PRECLS,TT_LSC,QT_LSC)
195 CALL LSCOND (PSG,dpFac,QG1,QSAT,
196 O PRECLS(1,myThid),TT_LSC(1,1,myThid),
197 O QT_LSC(1,1,myThid),
198 I kGround,bi,bj,myThid)
199
200 IF ( aim_energPrecip ) THEN
201 C 2.3 Snow Precipitation (update TT_CNV & TT_LSC)
202 CALL SNOW_PRECIP (
203 I PSG, dpFac, SE, ICLTOP,
204 I PRECNV(1,myThid), QT_CNV(1,1,myThid),
205 I PRECLS(1,myThid), QT_LSC(1,1,myThid),
206 U TT_CNV(1,1,myThid), TT_LSC(1,1,myThid),
207 O EnPrec(1,myThid),
208 I kGround,bi,bj,myThid)
209 ELSE
210 DO J=1,NGP
211 EnPrec(J,myThid) = 0. _d 0
212 ENDDO
213 ENDIF
214
215 C-- 3. Radiation (shortwave and longwave) and surface fluxes
216
217 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
218 C --> from FORDATE (in SPEEDY) :
219
220 C 3.0 Compute Incomming shortwave rad. (from FORDATE in SPEEDY)
221
222 c_FM CALL SOL_OZ (SOLC,TYEAR)
223 CALL SOL_OZ (SOLC,tYear, snLat(1,myThid), csLat(1,myThid),
224 O FSOL, OZONE, OZUPP, ZENIT, STRATZ,
225 I bi,bj,myThid)
226
227 C <-- from FORDATE (in SPEEDY).
228 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
229
230 C 3.1 Compute shortwave tendencies and initialize lw transmissivity
231
232 C Set LW absorption in CO2 band
233 IF ( aim_select_pCO2.EQ.1 .OR. aim_select_pCO2.EQ.3 ) THEN
234 absLW_CO2 = ABLCO2
235 & + aim_abs_pCO2*LOG( aim_pCO2/aim_ref_pCO2 )
236 ELSE
237 absLW_CO2 = ABLCO2
238 ENDIF
239
240 C The sw radiation may be called at selected time steps
241 LRADSW = .TRUE.
242
243 IF (LRADSW) THEN
244
245 c_FM CALL RADSW (PSG,QG1,RH,ALB1,
246 c_FM & ICLTOP,CLOUDC,TSR,SSR,TT_RSW)
247 ICLTOP(1) = 1
248 CALL RADSW (PSG,dpFac,QG1,RH(1,1,myThid),ALB1(1,0,myThid),
249 I FSOL, OZONE, OZUPP, ZENIT, STRATZ,
250 O TAU2, STRATC,
251 O ICLTOP,CLOUDC(1,myThid),
252 O TSR(1,myThid),SSR(1,0,myThid),
253 O UPSWG,TT_RSW(1,1,myThid),
254 I absLW_CO2, kGround, bi, bj, myThid )
255
256 DO J=1,NGP
257 CLTOP(J,myThid)=SIGH(ICLTOP(J)-1)*PSG_1(J)
258 ENDDO
259
260 DO K=1,NLEV
261 DO J=1,NGP
262 TT_RSW(J,K,myThid)=TT_RSW(J,K,myThid)*RPS_1*GRDSCP(K)
263 ENDDO
264 ENDDO
265
266 #ifdef ALLOW_DIAGNOSTICS
267 IF ( usePkgDiag ) THEN
268 CALL DIAGNOSTICS_FILL( UPSWG,
269 & 'UPSWG ', 1, 1 , 3,bi,bj, myThid )
270 ENDIF
271 #endif
272
273 ENDIF
274
275 C 3.2 Compute downward longwave fluxes
276
277 c_FM CALL RADLW (-1,TG1,TS,ST4S,
278 c_FM & OLR,SLR,TT_RLW)
279 CALL RADLW (-1,TG1,TS(1,myThid),ST4S,
280 & OZUPP, STRATC, TAU2, FLUX, ST4A,
281 O OLR(1,myThid),SLR(1,0,myThid),TT_RLW(1,1,myThid),
282 I kGround,bi,bj,myThid)
283
284 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
285 C 3.3. Compute surface fluxes and land skin temperature
286
287 c_FM CALL SUFLUX (PSG,UG1,VG1,TG1,QG1,RH,PHIG1,
288 c_FM & PHIS0,FMASK1,STL1,SST1,SOILW1,SSR,SLR,
289 c_FM & USTR,VSTR,SHF,EVAP,ST4S,
290 c_FM & TS,TSKIN,U0,V0,T0,Q0)
291 CALL SUFLUX_PREP(
292 I PSG, TG1, QG1, RH(1,1,myThid), SE, VsurfSq,
293 I WVSurf(1,myThid),csLat(1,myThid),fOrogr(1,myThid),
294 I FMASK1(1,1,myThid),STL1(1,myThid),SST1(1,myThid),
295 I sti1(1,myThid), SSR(1,0,myThid),
296 O SPEED0(1,myThid),DRAG(1,0,myThid),DENVV,
297 O dTskin,T1s,T0(1,myThid),Q0(1,myThid),
298 I kGround,bi,bj,myThid)
299
300 CALL SUFLUX_LAND (
301 I PSG, FMASK1(1,1,myThid), EMISFC,
302 I STL1(1,myThid), dTskin,
303 I SOILW1(1,myThid), SSR(1,1,myThid), SLR(1,0,myThid),
304 I T1s, T0(1,myThid), Q0(1,myThid), DENVV,
305 O SHF(1,1,myThid), EVAP(1,1,myThid), SLR(1,1,myThid),
306 O Shf0, dShf, Evp0, dEvp, Slr0, dSlr, sFlx,
307 O TS(1,myThid), TSKIN(1,myThid),
308 I bi,bj,myThid)
309 #ifdef ALLOW_LAND
310 CALL AIM_LAND_IMPL(
311 I FMASK1(1,1,myThid), dTskin,
312 I Shf0, dShf, Evp0, dEvp, Slr0, dSlr,
313 U sFlx, STL1(1,myThid),
314 U SHF(1,1,myThid), EVAP(1,1,myThid), SLR(1,1,myThid),
315 O dTsurf(1,1,myThid),
316 I bi, bj, myTime, myIter, myThid)
317 #endif /* ALLOW_LAND */
318
319 CALL SUFLUX_OCEAN(
320 I PSG, FMASK1(1,2,myThid),
321 I SST1(1,myThid),
322 I SSR(1,2,myThid), SLR(1,0,myThid),
323 O T1s, T0(1,myThid), Q0(1,myThid), DENVV,
324 O SHF(1,2,myThid), EVAP(1,2,myThid), SLR(1,2,myThid),
325 I bi,bj,myThid)
326
327 IF ( aim_splitSIOsFx ) THEN
328 CALL SUFLUX_SICE (
329 I PSG, FMASK1(1,3,myThid), EMISFC,
330 I STI1(1,myThid), dTskin,
331 I SSR(1,3,myThid), SLR(1,0,myThid),
332 I T1s, T0(1,myThid), Q0(1,myThid), DENVV,
333 O SHF(1,3,myThid), EVAP(1,3,myThid), SLR(1,3,myThid),
334 O Shf0, dShf, Evp0, dEvp, Slr0, dSlr, sFlx,
335 O TS(1,myThid), TSKIN(1,myThid),
336 I bi,bj,myThid)
337 #ifdef ALLOW_THSICE
338 CALL AIM_SICE_IMPL(
339 I FMASK1(1,3,myThid), SSR(1,3,myThid), sFlx,
340 I Shf0, dShf, Evp0, dEvp, Slr0, dSlr,
341 U STI1(1,myThid),
342 U SHF(1,3,myThid), EVAP(1,3,myThid), SLR(1,3,myThid),
343 O dTsurf(1,3,myThid),
344 I bi, bj, myTime, myIter, myThid)
345 #endif /* ALLOW_THSICE */
346 ELSE
347 DO J=1,NGP
348 SHF (J,3,myThid) = 0. _d 0
349 EVAP(J,3,myThid) = 0. _d 0
350 SLR (J,3,myThid) = 0. _d 0
351 ENDDO
352 ENDIF
353
354 CALL SUFLUX_POST(
355 I FMASK1(1,1,myThid), EMISFC,
356 I STL1(1,myThid), SST1(1,myThid), sti1(1,myThid),
357 I dTskin, SLR(1,0,myThid),
358 I T0(1,myThid), Q0(1,myThid), DENVV,
359 U DRAG(1,0,myThid), SHF(1,0,myThid),
360 U EVAP(1,0,myThid), SLR(1,1,myThid),
361 O ST4S, TS(1,myThid), TSKIN(1,myThid),
362 I bi,bj,myThid)
363
364 #ifdef ALLOW_DIAGNOSTICS
365 IF ( usePkgDiag ) THEN
366 CALL DIAGNOSTICS_FILL( SLR(1,0,myThid),
367 & 'DWNLWG ', 1, 1 , 3,bi,bj, myThid )
368 ENDIF
369 #endif
370 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
371
372 C 3.4 Compute upward longwave fluxes, convert them to tendencies
373 C and add shortwave tendencies
374
375 c_FM CALL RADLW (1,TG1,TS,ST4S,
376 c_FM & OLR,SLR,TT_RLW)
377 CALL RADLW (1,TG1,TS(1,myThid),ST4S,
378 & OZUPP, STRATC, TAU2, FLUX, ST4A,
379 O OLR(1,myThid),SLR(1,0,myThid),TT_RLW(1,1,myThid),
380 I kGround,bi,bj,myThid)
381
382 DO K=1,NLEV
383 DO J=1,NGP
384 TT_RLW(J,K,myThid)=TT_RLW(J,K,myThid)*RPS_1*GRDSCP(K)
385 c_FM TTEND (J,K)=TTEND(J,K)+TT_RSW(J,K)+TT_RLW(J,K)
386 ENDDO
387 ENDDO
388
389 #ifdef ALLOW_CLR_SKY_DIAG
390 C 3.5 Compute clear-sky radiation (for diagnostics only)
391 IF ( aim_clrSkyDiag ) THEN
392
393 C 3.5.1 Compute shortwave tendencies
394 dummyI(1) = -1
395 CALL RADSW (PSG,dpFac,QG1,RH(1,1,myThid),ALB1(1,0,myThid),
396 I FSOL, OZONE, OZUPP, ZENIT, STRATZ,
397 O TAU2, STRATC,
398 O dummyI, dummyR,
399 O TSWclr(1,myThid), SSWclr(1,myThid), UPSWG, TT_SWclr(1,1,myThid),
400 I absLW_CO2, kGround, bi, bj, myThid )
401
402 #ifdef ALLOW_DIAGNOSTICS
403 IF ( usePkgDiag ) THEN
404 CALL DIAGNOSTICS_FILL( UPSWG,
405 & 'UPSWGclr', 1, 1 , 3,bi,bj, myThid )
406 ENDIF
407 #endif
408
409 C 3.5.2 Compute downward longwave fluxes
410
411 CALL RADLW (-1,TG1,TS(1,myThid),ST4S,
412 & OZUPP, STRATC, TAU2, FLUX, ST4A,
413 O OLWclr(1,myThid), SLWclr(1,myThid), TT_LWclr(1,1,myThid),
414 I kGround,bi,bj,myThid)
415
416 C 3.5.3 Compute upward longwave fluxes, convert them to tendencies
417
418 CALL RADLW (1,TG1,TS(1,myThid),ST4S,
419 & OZUPP, STRATC, TAU2, FLUX, ST4A,
420 O OLWclr(1,myThid), SLWclr(1,myThid), TT_LWclr(1,1,myThid),
421 I kGround,bi,bj,myThid)
422
423 DO K=1,NLEV
424 DO J=1,NGP
425 TT_SWclr(J,K,myThid)=TT_SWclr(J,K,myThid)*RPS_1*GRDSCP(K)
426 TT_LWclr(J,K,myThid)=TT_LWclr(J,K,myThid)*RPS_1*GRDSCP(K)
427 ENDDO
428 ENDDO
429
430 ENDIF
431 #endif /* ALLOW_CLR_SKY_DIAG */
432
433 C-- 4. PBL interactions with lower troposphere
434
435 C 4.1 Vertical diffusion and shallow convection
436
437 c_FM CALL VDIFSC (UG1,VG1,SE,RH,QG1,QSAT,PHIG1,
438 c_FM & UT_PBL,VT_PBL,TT_PBL,QT_PBL)
439 CALL VDIFSC (dpFac, SE, RH(1,1,myThid), QG1, QSAT,
440 O TT_PBL(1,1,myThid),QT_PBL(1,1,myThid),
441 I kGround,bi,bj,myThid)
442
443 C 4.2 Add tendencies due to surface fluxes
444
445 DO J=1,NGP
446 c_FM UT_PBL(J,NLEV)=UT_PBL(J,NLEV)+USTR(J,3)*RPS(J)*GRDSIG(NLEV)
447 c_FM VT_PBL(J,NLEV)=VT_PBL(J,NLEV)+VSTR(J,3)*RPS(J)*GRDSIG(NLEV)
448 c_FM TT_PBL(J,NLEV)=TT_PBL(J,NLEV)+ SHF(J,3)*RPS(J)*GRDSCP(NLEV)
449 c_FM QT_PBL(J,NLEV)=QT_PBL(J,NLEV)+EVAP(J,3)*RPS(J)*GRDSIG(NLEV)
450 K = kGround(J)
451 IF ( K.GT.0 ) THEN
452 TT_PBL(J,K,myThid) = TT_PBL(J,K,myThid)
453 & + SHF(J,0,myThid) *RPS_1*GRDSCP(K)
454 QT_PBL(J,K,myThid) = QT_PBL(J,K,myThid)
455 & + EVAP(J,0,myThid)*RPS_1*GRDSIG(K)
456 ENDIF
457 ENDDO
458
459 c_FM DO K=1,NLEV
460 c_FM DO J=1,NGP
461 c_FM UTEND(J,K)=UTEND(J,K)+UT_PBL(J,K)
462 c_FM VTEND(J,K)=VTEND(J,K)+VT_PBL(J,K)
463 c_FM TTEND(J,K)=TTEND(J,K)+TT_PBL(J,K)
464 c_FM QTEND(J,K)=QTEND(J,K)+QT_PBL(J,K)
465 c_FM ENDDO
466 c_FM ENDDO
467
468 #endif /* ALLOW_AIM */
469
470 RETURN
471 END

  ViewVC Help
Powered by ViewVC 1.1.22