/[MITgcm]/MITgcm_contrib/AITCZ/code/mitphys_do_atmos_physics.F
ViewVC logotype

Contents of /MITgcm_contrib/AITCZ/code/mitphys_do_atmos_physics.F

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


Revision 1.1 - (show annotations) (download)
Wed Aug 20 15:24:59 2003 UTC (21 years, 11 months ago) by czaja
Branch: MAIN
CVS Tags: HEAD
Initial creation of Arnaud's simple coupled simulation.

1 C $Header: /u/gcmpack/models/MITgcmUV/verification/aim.5l_Equatorial_Channel/code/aim_do_atmos_physics.F,v 1.4.2.1 2001/04/30 23:26:33 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5 #undef OLD_AIM_GRIG_MAPPING
6
7 SUBROUTINE MITPHYS_DO_ATMOS_PHYSICS( phi_hyd, currentTime, myThid)
8 C /==================================================================\
9 C | S/R MITPHYS_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 Physics package
23 #include "atparam.h"
24 #include "atparam1.h"
25 #include "MITPHYS_PARAMS.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 #include "thermo.h"
32 #include "com_mitphysvar.h"
33 #include "com_forcing1.h"
34 #include "Lev_def.h"
35 C MITgcm
36 #include "EEPARAMS.h"
37 #include "PARAMS.h"
38 #include "DYNVARS.h"
39 #include "GRID.h"
40 #include "SURFACE.h"
41
42 C -------------- Routine arguments -----------------------------------
43 _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
44 _RL currentTime
45 INTEGER myThid
46
47 C omp modif: shapiro filter on the convective mass flux:
48 C _RL cbmf_tmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
49 _RL dist
50
51
52 #ifdef ALLOW_MITPHYS
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 _RL 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 _RL 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 _RL hInitial(Nr)
93 _RL hInitialW(Nr)
94 _RL pSurfs(Nr)
95 c 5 levels DATA pSurfs / 950.D2,775.D2, 500.D2, 250.D2, 75.D2 /
96 DATA pSurfs / 1000.0D2, 975.0D2, 950.0D2, 925.0D2, 900.0D2,
97 : 875.0D2, 850.0D2, 825.0D2, 800.0D2, 775.0D2,
98 : 750.0D2, 725.0D2, 700.0D2, 675.0D2, 650.0D2,
99 : 625.0D2, 600.0D2, 575.0D2, 550.0D2, 525.0D2,
100 : 500.0D2, 475.0D2, 450.0D2, 425.0D2, 400.0D2,
101 : 375.0D2, 350.0D2, 325.0D2, 300.0D2, 275.0D2,
102 : 250.0D2, 225.0D2, 200.0D2, 175.0D2, 150.0D2,
103 : 125.0D2, 100.0D2, 75.0D2, 50.0D2, 25.0D2 /
104 SAVE pSurfs
105 _RL Soilqmax
106 _RL pGround
107 INTEGER bi, bj
108 INTEGER Katm
109
110 INTEGER NEXT_CALL
111 DATA NEXT_CALL /0/
112 SAVE NEXT_CALL
113
114
115 _RL LAT_CYCLE
116
117 C LOGICAL DO_OMP_SHAP
118 C
119
120 CC (acz) check for SST
121 CC write(6,*) 'beginning1 of do_atmos_phys: SST_REF=',SST_REF(1)
122
123 bi = 1
124 bj = 1
125
126 IF (CurrentTime .LT. 10.d0) THEN
127 DO J=1,sNy
128 DO I=1,sNx
129 I2 = (sNx)*(J-1)+I
130 SST1 (I2)= SST_REF(I2)
131 STL1 (I2)= STL_REF(I2)
132 CBMF_DYN(I,J,bi,bj) = CBMFG1(I2)
133 SST_DYN(I,J,bi,bj) = SST_REF(I2)
134 STL_DYN(I,J,bi,bj) = STL_REF(I2)
135 END DO
136 END DO
137 END IF
138
139 CC (acz) check for SST
140 CC write(6,*) 'beginning2 of do_atmos_phys: SST1=',SST1(1)
141 CC write(6,*) 'beginning2 of do_atmos_phys: SST_DYN=',SST_DYN(1,1,bi,bj)
142
143 IF (NEXT_CALL .EQ. 0 ) THEN
144
145 pGround = 1012.5D2
146
147 C Assume only one tile per proc. for now
148 bi = 1
149 bj = 1
150 IG0 = myXGlobalLo
151 JG0 = myYGlobalLo
152
153 C
154 C Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
155 C Internal index mapping is linear in X and Y with a second
156 C dimension for the vertical.
157
158 C Adjustment for heave due to mean heating/cooling
159 C ( I don't think the old formula was strictly "correct" for orography
160 C but I have implemented it as was for now. As implemented
161 C the mean heave of the bottom (K=Nr) level is calculated rather than
162 C the mean heave of the base of the atmosphere. )
163
164 c -- sb: remove this with the MIT physics:
165 c sb c_jmc: Because AIM physics LSC is not applied in the stratosphere (top level),
166 c sb c ==> move water wapor from the stratos to the surface level.
167 c sb DO J = 1-Oly, sNy+Oly
168 c sb DO I = 1-Olx, sNx+Olx
169 c sb c k = k_surf(i,j,bi,bj)
170 c sb c salt(I,J,k,bi,bj) = salt(I,J,k,bi,bj)
171 c sb c & + maskC(i,j,Nr,bi,bj)*salt(I,J,Nr,bi,bj)*drF(Nr)*recip_drF(k)
172 c sb salt(I,J,Nr,bi,bj) = 0.
173 c sb ENDDO
174 c sb ENDDO
175 c sb --
176
177 C Note the mapping here is only valid for one tile per proc.
178 DO J = 1, sNy
179 DO I = 1, sNx
180 I2 = (sNx)*(J-1)+I
181 C omp: done in mitphys_do_inphys now
182 C latitude and logitude, assuming a spherical coordinate system
183 C Must be modified for cartesian coordinates
184 LAT_G(I2) = yC(I,J,bi,bj)
185 LON_G(I2) = xC(I,J,bi,bj)
186 C omp: new addition for restart files.
187 CBMFG1(I2) = CBMF_DYN (I,J,bi, bj)
188 SST1(I2) = SST_DYN (I,J,bi, bj)
189 STL1(I2) = STL_DYN (I,J,bi, bj)
190
191 CC (acz) check for SST
192 CC IF (I2 .EQ. 1) THEN
193 CC write(6,*) 'SST1 set to SST_DYN at the begin. of do_atmos_phys: SST1,SSTDYN', SST1(I2),SST_DYN (1,1,bi, bj)
194 CC END IF
195
196
197 DO K = 1, Nr
198 Katm = K
199
200 C UG1(I2,Katm) = 0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
201 C VG1(I2,Katm) = 0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))
202
203 C UG1(I2,Katm) = uVel(I,J,K,bi,bj)
204 C VG1(I2,Katm) = vVel(I,J,K,bi,bj)
205
206 UGW(I2,Katm) = uVel(I,J,K,bi,bj)
207 UGE(I2,Katm) = uVel(I+1,J,K,bi,bj)
208 VGS(I2,Katm) = vVel(I,J,K,bi,bj)
209 VGN(I2,Katm) = vVel(I,J+1,K,bi,bj)
210
211 WG1(I2, Katm) = wVel (I,J,K,bi, bj)
212
213 C Physics works with temperature - not potential temp.
214 TG1(I2,Katm) = theta(I,J,K,bi,bj)/
215 & ((pGround/pSurfs(K))**(RD/CPD))
216 c_jmc QG1(I2,Katm) = salt(I,J,K,bi,bj)
217 QG1(I2,Katm) = MAX(salt(I,J,K,bi,bj), 0. _d 0)
218 PHIG1(I2,Katm) = phi_hyd(I,J,K,bi,bj) + etaN(I,J,bi,bj) ! not used in MIT physics
219
220 ENDDO
221
222 ENDDO
223 ENDDO
224
225 C
226 DO J=1,sNy
227 DO I=1,sNx
228 I2 = (sNx)*(J-1)+I
229 PSG1(I2) = pGround
230 ENDDO
231 ENDDO
232 C write(*,*) 'CBMF_DYN before:', CBMF_DYN
233 C write(*,*) 'CBMFGE before:', CBMFG1
234
235 C
236 C
237 C Physics package needs to know time of year as a fraction
238
239 C not used anymore -omp
240
241 tYear = currentTime/(86400.*360.) -
242 & FLOAT(INT(currentTime/(86400.*360.)))
243 C_EqCh: Fix solar insolation with Sun directly overhead on Equator
244 tYear = 0.25
245 C
246 C Load external data needed by physics package
247 C 1. Albedo
248 C 2. Surface temperatures
249 C 3. Soil moisture
250 C 4. Snow depth - assume no snow for now
251 C 5. Sea ice - assume no sea ice for now
252 C 6. Land sea mask - infer from exact zeros in soil moisture dataset
253 C 7. Surface geopotential - to be done when orography is in
254 C dynamical kernel. Assume 0. for now.
255 mnthIndex = INT(tYear*12.)+1
256 C IF ( mnthIndex .NE. prevMnthIndex .OR.
257 C & FirstCall ) THEN
258 C prevMnthIndex = mnthIndex
259
260 C 1. albedo (not used)
261
262 DO J=1,sNy
263 DO I=1,sNx
264 I2 = (sNx)*(J-1)+I
265 alb0(I2) = 0.
266
267 ENDDO
268 ENDDO
269
270 C modif omp: introduce an annual cycle in the prescribed SST.
271 C modif acz: introduce an annual cycle in the prescribed STL.
272
273 IF (ANNUAL_CYCLE) THEN
274
275 C (acz) LAT_CYCLE = 1 at t=0 to be consistent with do_inphys.F
276 C This means that the integration starts in late summer
277
278 LAT_CYCLE = 0.5 + 0.5 * cos(2 * pi *
279 & (currentTime / 24./ 3600.) / 360.)
280
281
282 C 2. SST / STL
283 C modif omp: the background SST is defined thourgh the MITPHYS namelist
284 C modif acz: introduce analytical expression for STL
285
286 C omp: for spherical geometry
287 IF (UsingSphericalPolarGrid) THEN
288
289 DO J=1,sNy
290 DO I=1,sNx
291 I2 = (sNx)*(J-1)+I
292
293 SST_REF(I2) = SST_BACK
294 STL_REF(I2) = SST_BACK
295
296 C Equator to pole SST gradient
297
298 SST_REF(I2) = SST_REF(I2) - DELT_EQ_PO *
299 & sin( LAT_G(I2) * pi / 180. ) **2
300
301
302 C Equator to northern boundary SST gradient - work only for a single tile
303
304 SST_REF(I2) = SST_REF(I2) - DELT_EQ_BND *
305 & sin( LAT_G(I2) / abs(Phimin)) **2
306
307 C Lindzen and Hou type gradient:
308
309 SST_REF(I2) = SST_REF(I2) - DELT_LH * (
310 & 3 * ( sin ( PI * LAT_G(I2) / 180.d0 )
311 & - sin(PI * LAT_CYCLE * LAT_LH / 180.d0)
312 & ) ** 2)
313
314
315 CC (acz) SST used for TOSA1:
316 CC Zonally symmetric SST perturbation
317
318
319 dist = SQRT ( (LAT_G(I2) - LAT_PERT*LAT_CYCLE) ** 2)
320
321 if (dist .LE. DEL_LAT) then
322 SST_REF(I2) = SST_REF(I2) + DELT_PERT
323 & * cos ( pi * dist / 2./DEL_LAT ) **2
324 end if
325
326 CC (acz) STL used for TOSA1:
327 CC
328 CC 1 - Warm ring centered at the equator with same
329 CC meridional extent (DEL_LAT) than the above (TOSA1) SST
330 CC The profile is further time-dependent:
331 C strongest at equinox, zero at solstice
332
333 dist = ABS ( LAT_G(I2) )
334
335 if (dist .LE. DEL_LAT) then
336 STL_REF(I2) = STL_REF(I2) + DELT_STL_EQU
337 & * cos ( pi * dist / 2./DEL_LAT ) **2
338 & * cos ( 2 * pi * (currentTime / 24./ 3600.)
339 & / 360. + 0.5 * pi ) **2
340 end if
341
342 C 2- Opposite sign variations north and south of equator
343 C North is warm at t=0
344
345 dist = ABS ( LAT_G(I2) )
346
347 if (dist .LE. DEL_LAT_STL) then
348 STL_REF(I2) = STL_REF(I2) + DELT_STL_HEM
349 & * sin ( pi * LAT_G(I2) / DEL_LAT_STL )
350 & * cos ( 2 * pi * (currentTime / 24./ 3600.) / 360.)
351 end if
352
353 ENDDO
354 ENDDO
355
356 ELSE
357
358 C Cartesian Grid
359 DO J=1,sNy
360 DO I=1,sNx
361 I2 = (sNx)*(J-1)+I
362 SST_REF(I2) = SST_BACK
363 SST_REF(I2) = SST_REF(I2) - DELT_EQ_BND *
364 & sin( 0.5 * PI * yC(I,J,bi,bj) / abs(phiMin) ) **2
365
366 END DO
367 END DO
368
369 END IF
370
371 END IF
372 CC end for the ANNUAL_CYCLE end if loop
373
374 IF (PRESC_SST) THEN
375 DO I2=1,NGP
376 SST1(I2) = SST_REF(I2)
377 END DO
378 END IF
379
380 IF (PRESC_STL) THEN
381 DO I2=1,NGP
382 STL1(I2) = STL_REF(I2)
383 END DO
384 END IF
385
386 CC (acz) check for SST
387 CC write(6,*) 'end of do_atmos_physic: SST1=',SST1(1)
388
389
390 C 3. Soil moisture (not used)
391
392 DO J=1,sNy
393 DO I=1,sNx
394 I2 = (sNx)*(J-1)+I
395 soilq1(I2) = 0.
396 ENDDO
397 ENDDO
398
399 Soilqmax=20.
400
401 C ENDIF
402 C
403
404 C 4. Snow depth (not used)
405 IF ( FirstCall ) THEN
406 C Set snow depth, sea ice to zero for now
407 C Land-sea mask ( figure this out from where
408 C soil moisture is exactly zero ).
409 CC
410 CC(acz) VERY RISKY lines !!!!!!!!!
411 CC DO J=1,sNy
412 CC DO I=1,sNx
413 CC I2 = (sNx)*(J-1)+I
414 CC fMask1(I2) = 1.
415 CC IF ( soilq1(I2) .EQ. 0. ) fMask1(I2) = 0.
416 CC oice1(I2) = 0.
417 CC snow1(I2) = 0.
418 CC ENDDO
419 CC ENDDO
420
421 ENDIF
422 C
423 C Addition may 15 . Reset humidity to 0. if negative
424 C ---------------------------------------------------
425 Caja DO K=1,Nr
426 Caja DO J=1-OLy,sNy+OLy
427 Caja DO I=1-Olx,sNx+OLx
428 Caja IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
429 Caja salt(i,j,k,bi,bj) = 0.
430 Caja ENDIF
431 Caja ENDDO
432 Caja ENDDO
433 Caja ENDDO
434
435 cccc CALL PDRIVER( tYear )
436 CCC (acz) also known as mitphy_driver.F !!!
437 CCC
438 write(*,*) 'currentTime: ',currentTime
439 CALL MPDRIVER( currentTime, float(PHYS_CALL) * deltaT ) ! Call the MIT physics
440 C
441 #ifdef ALLOW_TIMEAVE
442 C Calculate diagnostics for MITPHYS
443 CALL MITPHYS_CALC_DIAGS( bi, bj, currentTime, myThid )
444 #endif /* ALLOW_TIMEAVE */
445 C
446 FirstCall = .FALSE.
447
448 C
449 C The dry adiabatic adjustment done in the convection scheme may
450 C have affected the temperature and humidity profiles:
451 C
452 DO K = 1, Nr
453 DO J = 1, sNy
454 c! DO J = 1, sNy-1 ! because last row of latitude not meaningful here
455 DO I = 1, sNx
456 I2 = (sNx)*(J-1)+I
457 Katm = K
458 theta(I,J,K,bi,bj) = TG1(I2,Katm)
459 : * ( (pGround/pSurfs(K))**(RD/CPD) )
460 salt(I,J,K,bi,bj) = QG1(I2,Katm)
461 C omp: dry adjustment on velocity is now done as a time-tendency term in external_forcing.
462 C uVel(I,J,K,bi,bj) = UGW(I2,Katm)
463 C vVel(I,J,K,bi,bj) = VGS(I2,Katm)
464 ENDDO
465 ENDDO
466 END DO
467
468
469 DO J = 1, sNy
470 DO I = 1, sNx
471 I2 = (sNx)*(J-1)+I
472 CBMF_DYN(I,J,bi,bj) = CBMFG1(I2)
473 SST_DYN(I,J,bi,bj) = SST1(I2)
474 STL_DYN(I,J,bi,bj) = STL1(I2)
475 END DO
476 END DO
477
478 C -- Following section is neccesary since we are adjusting theta/salt
479 C directly and the overlaps need to be updated accordingly.
480 _EXCH_XYZ_R8(theta,myThid)
481 _EXCH_XYZ_R8(salt,myThid)
482 _EXCH_XY_R8(CBMF_DYN,myThid)
483 _EXCH_XY_R8(SST_DYN,myThid)
484 _EXCH_XY_R8(STL_DYN,myThid)
485
486
487 C omp: shapiro-filtering the convective mass flux:
488 C -still assuming single processor:
489 C omp: dropped... to be implemented again.
490
491 c$$$
492 c$$$ DO_OMP_SHAP = .FALSE.
493 c$$$ IF (DO_OMP_SHAP)THEN
494 c$$$ DO J = 1, sNy
495 c$$$ DO I = 1, sNx
496 c$$$ I2 = I+(J-1)*sNx
497 c$$$ cbmf_tmp(i,j,bi,bj) = CBMFG1(I2)
498 c$$$ end do
499 c$$$ end do
500 c$$$
501 c$$$ CALL OMP_SHAP_FILTER(cbmf_tmp,currentTime, myThid)
502 c$$$
503 c$$$ DO J = 1, sNy
504 c$$$ DO I = 1, sNx
505 c$$$ I2 = I+(J-1)*sNx
506 c$$$ CBMFG1(I2) = CBMF_tmp(i,j,bi,bj)
507 c$$$ end do
508 c$$$ end do
509 c$$$ END IF
510
511 C
512 C
513 C For velocity tendencies on C grid need to interpolate
514 C to do this in multi-processor mode we copy U and V tendencies
515 C into dynamics conforming array and then exchange.
516 C ** NOTE - Exchange at this point is not compatible with
517 C multiple-tiles per compute thread/process.
518 C However, AIM code itself is not compatible with
519 C this as its global data assumes only one "tile".
520 c sb CALL AIM_UTVT2DYN( bi, bj, myThid )
521 c sb _EXCH_XYZ_R8( AIM_UT, myThid )
522 c sb _EXCH_XYZ_R8( AIM_VT, myThid )
523
524 c sb: remove this for the moment (no drag computed by MITPHYS):
525 c sb DO J = 1, sNy
526 c sb DO I = 1, sNx
527 c sb I2 = I+(J-1)*sNx
528 c sb aim_drag(I,J,bi,bj) = DRAG(I2)
529 c sb ENDDO
530 c sb ENDDO
531 c sb _EXCH_XY_R8( aim_drag, myThid )
532 C
533 NEXT_CALL = PHYS_CALL
534 C write(*,*) 'CBMF_DYN after:', CBMF_DYN
535 c write(*,*) 'CBMFG1 after:', CBMFG1
536 C write(*,*) 'LAT :', LAT_G
537 C write(*,*) 'LON :', LON_G
538 END IF
539 C(acz) end if for the NEXT_CALL if
540 NEXT_CALL = NEXT_CALL -1
541
542 #endif /* ALLOW_MITPHYS */
543
544 RETURN
545 END
546
547
548

  ViewVC Help
Powered by ViewVC 1.1.22