source: trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F

Last change on this file was 2874, checked in by jleconte, 23 months ago

Changed nomenclature for planet variables from m_ to p_ in interface V4
File size: 27.2 KB
Line 
1!*******************************************************************************
2! PURPOSE: LMD_driver is the WRF mediation layer routine that provides the
3! interface to LMD physics packages in the WRF model layer. For those familiar
4! with the LMD GCM, the aim of this driver is to do part of the job of the
5! calfis.F routine.
6!*******************************************************************************
7! AUTHOR: A. Spiga - January 2007
8!*******************************************************************************
9! UPDATES:  - included all final updates for the paper - March 2008
10!           - general cleaning of code and comments - October 2008
11!           - additions for idealized cases - January 2009
12!           - additions for new soil model in physics - January 2010
13!           - unified module_lmd_driver: old, new phys and LES - February 2011
14!           - a new paradigm to prepare nested simulations - April 2014
15!           - adapted to new interface philosophy & other planets - August 2016
16!           - adapated to WRFV4 - JL22
17!*******************************************************************************
18MODULE module_lmd_driver
19CONTAINS
20    SUBROUTINE lmd_driver(id,max_dom,DT,ITIMESTEP,XLAT,XLONG, &
21        IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME,KMS,KME, &
22        i_start,i_end,j_start,j_end,kts,kte,num_tiles, &
23        DX,DY, &
24        MSFT,MSFU,MSFV, &
25        GMT,JULYR,JULDAY, &
26        P8W,DZ8W,T8W,Z,HT,MUT,DNW, &
27        U,V,TH,T,P,EXNER,RHO, &
28        P_HYD, P_HYD_W, &
29        PTOP, &
30        RADT, &
31        TSK,PSFC, &
32        RTHPLATEN,RUPLATEN,RVPLATEN, &
33        num_3d_s,SCALAR, &
34        num_3d_m,moist, &
35        TRACER_MODE, &
36        planet_type, &
37        P_ALBEDO,P_TI,P_CO2ICE,P_EMISS, &
38        P_H2OICE,P_TSOIL,P_Q2,P_TSURF, &
39        P_FLUXRAD,P_WSTAR,P_ISOIL,P_DSOIL,&
40        P_Z0, CST_Z0, P_GW, &
41        NUM_SOIL_LAYERS, &
42        CST_AL, CST_TI, &
43        isfflx, diff_opt, km_opt, &
44        HR_SW,HR_LW,HR_DYN,DT_RAD,&
45        CLOUDFRAC,TOTCLOUDFRAC,RH, &
46        DQICE,DQVAP,REEVAP,SURFRAIN,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
47        FLUXSURF_LW,FLXGRD,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF,LATENT_HF,&
48        HFMAX,ZMAX,&
49        USTM,HFX,&
50        SLPX,SLPY,RESTART)
51! NB: module_lmd_driver_output1.inc : output arguments generated from Registry
52
53
54
55
56
57!==================================================================
58! USES
59!==================================================================
60   USE module_state_description
61   USE module_model_constants
62   USE module_wrf_error
63   !!!!!!!! interface modules
64   USE variables_mod !! to set variables
65   USE update_inputs_physiq_mod !! to set inputs for physiq
66   USE update_outputs_physiq_mod !! to get outputs from physiq
67   USE iniphysiq_mod !! to get iniphysiq subroutine
68   USE callphysiq_mod !! to call the LMD physics
69   !!!!!!!! interface modules
70
71!==================================================================
72  IMPLICIT NONE
73!==================================================================
74
75!==================================================================
76! VARIABLES
77!==================================================================
78
79CHARACTER(len=10),INTENT(IN) :: planet_type
80! WRF Dimensions
81INTEGER, INTENT(IN   )    ::   &
82      ids,ide,jds,jde,kds,kde, &
83      ims,ime,jms,jme,kms,kme, &
84            kts,kte,num_tiles, &
85              NUM_SOIL_LAYERS
86INTEGER, DIMENSION(num_tiles), INTENT(IN) ::  &
87      i_start,i_end,j_start,j_end   
88! Scalars
89INTEGER, INTENT(IN  ) :: JULDAY, itimestep, julyr,id,max_dom,TRACER_MODE
90INTEGER, INTENT(IN  ) :: isfflx,diff_opt,km_opt
91REAL, INTENT(IN  ) :: GMT,dt,dx,dy,RADT
92REAL, INTENT(IN  ) :: CST_AL, CST_TI
93REAL, INTENT(IN  ) :: PTOP
94! 2D arrays         
95REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: &
96     MSFT,MSFU,MSFV, &
97     XLAT,XLONG,HT,  &
98     MUT, & !total dry air column mass (in Pa)
99     P_ALBEDO,P_TI,P_EMISS, &
100     SLPX,SLPY, &
101     P_CO2ICE,P_H2OICE, &
102     P_TSURF, P_Z0, &
103     P_FLUXRAD,P_WSTAR, &
104     PSFC,TSK
105REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: &
106     HFMAX,ZMAX,&
107     USTM,HFX,TOTCLOUDFRAC,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
108     FLUXSURF_LW,FLXGRD,LATENT_HF,REEVAP,SURFRAIN
109REAL, DIMENSION(kms:kme), INTENT(IN ) :: DNW  ! del(eta), depth between full levels in eta variables.
110! 3D arrays
111REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
112     dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th,p_hyd,p_hyd_w
113REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
114     RTHPLATEN,RUPLATEN,RVPLATEN, &
115     HR_SW,HR_LW,HR_DYN,DT_RAD,&
116     CLOUDFRAC,RH,DQICE,DQVAP,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF
117REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: &
118     P_Q2
119REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(INOUT )  :: &
120     P_TSOIL,P_ISOIL, P_DSOIL
121REAL, INTENT(IN  ) :: CST_Z0
122REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: &
123     P_GW
124! 4D arrays
125INTEGER, INTENT(IN ) :: num_3d_s,num_3d_m
126REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_s ), INTENT(INOUT ) :: scalar
127REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_m ), INTENT(INOUT ) ::  moist
128! Logical
129LOGICAL, INTENT(IN ) :: restart
130
131!-------------------------------------------   
132! LOCAL  VARIABLES
133!-------------------------------------------
134   INTEGER ::    i,j,k,its,ite,jts,jte,ij,kk
135   INTEGER ::    subs,iii
136   
137
138   ! *** for tracer Mode 20 ***
139   REAL ::    tau_decay
140   ! *** ----------------------- ***
141
142   ! *** for LMD physics
143   ! ------> inputs:
144   INTEGER :: ngrid,nlayer,nq,nsoil
145   REAL :: MY 
146   REAL :: phisfi_val
147   LOGICAL :: lastcall
148   ! ----------
149
150   ! <------ outputs:
151   !     physical tendencies
152   ! ... intermediate arrays
153   REAL, DIMENSION(:), ALLOCATABLE  :: &
154     dz8w_prof,p8w_prof,p_prof,t_prof,t8w_prof, &
155     u_prof,v_prof,z_prof
156   REAL, DIMENSION(:,:), ALLOCATABLE :: q_prof
157
158   ! Additional control variables
159   INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe
160   REAL :: elaps
161   INTEGER :: test
162   REAL :: wappel_phys
163   LOGICAL, SAVE :: flag_first_restart
164   LOGICAL, SAVE :: firstcall = .true.
165   INTEGER, SAVE :: previous_id = 0
166   REAL, DIMENSION(:), ALLOCATABLE, SAVE :: dp_save
167   REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: du_save, dv_save, dt_save
168   REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq_save     
169
170!!!IDEALIZED IDEALIZED
171      REAL :: lat_input, lon_input, ls_input, lct_input
172      LOGICAL :: isles
173!!!IDEALIZED IDEALIZED
174
175   REAL :: tk1,tk2
176
177!==================================================================
178! CODE
179!==================================================================
180
181print *, '** ',planet_type,'** ENTER WRF-LMD DRIVER'
182
183!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184!! determine here if this is turbulence-resolving or not
185IF (JULYR .le. 8999) THEN
186    isles = .false.  ! "True" LES is not available in this version
187    IF (firstcall .EQV. .true.) PRINT *, '*** REAL-CASE SIMULATION ***'
188ELSE
189     IF (firstcall .EQV. .true.) PRINT *, '*** IDEALIZED SIMULATION ***'
190     IF ((diff_opt .eq. 2) .and. (km_opt .eq. 2)) THEN
191        IF (firstcall .EQV. .true.) THEN
192          PRINT *, '*** type: LES ***'
193          PRINT *, '*** diff_opt = 2 *** km_opt = 2'
194          PRINT *, '*** forcing is isfflx = ',isfflx
195        ENDIF
196        isles = .true.
197        !! SPECIAL LES
198     ELSE
199        IF (firstcall .EQV. .true.) THEN
200          PRINT *, '*** type: not LES ***'
201          PRINT *, '*** diff_opt = ',diff_opt
202          PRINT *, '*** km_opt = ',km_opt
203        ENDIF
204        isles = .false.
205        !! IDEALIZED, no LES
206        !! cependant, ne veut-on pas pouvoir
207        !! prescrire un flux en idealise ??
208     ENDIF
209ENDIF
210!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
211
212!-------------------------!
213! TWEAK on WRF DIMENSIONS !
214!-------------------------!
215its = i_start(1)   ! define here one big tile
216ite = i_end(num_tiles)
217jts = j_start(1)
218jte = j_end(num_tiles)
219!!
220IF (isles .eqv. .false.) THEN
221 relax=0
222 sponge_top=0               ! another value than 0 triggers instabilities 
223 IF (id .gt. 1) relax=2     ! fix to avoid noise in nesting simulations ; 1 >> too much noise ...
224ENDIF
225ips=its
226ipe=ite
227jps=jts
228jpe=jte
229IF (isles .eqv. .false.) THEN
230 IF (ips .eq. ids)   ips=its+relax !! IF tests necesary for parallel runs
231 IF (ipe .eq. ide-1) ipe=ite-relax
232 IF (jps .eq. jds)   jps=jts+relax
233 IF (jpe .eq. jde-1) jpe=jte-relax
234ENDIF
235kps=kts         !! start at surface
236IF (isles .eqv. .false.) THEN
237 kpe=kte-sponge_top
238ELSE
239  IF (firstcall .EQV. .true.) PRINT *, '*** IDEALIZED SIMULATION: LES *** kpe=kte'
240 kpe=kte !-sponge_top
241ENDIF
242
243!----------------------------!
244! DIMENSIONS FOR THE PHYSICS !
245!----------------------------!
246wappel_phys = RADT
247zdt_split = dt*wappel_phys            ! physical timestep (s)
248ngrid=(ipe-ips+1)*(jpe-jps+1)         ! size of the horizontal grid
249nlayer = kpe-kps+1                    ! number of vertical layers: nlayermx
250nsoil = NUM_SOIL_LAYERS               ! number of soil layers: nsoilmx
251CALL update_inputs_physiq_tracers(TRACER_MODE,nq)
252IF (firstcall .EQV. .true.) PRINT *,'** ',planet_type,'** TRACER MODE', TRACER_MODE
253
254! **** needed but hardcoded
255lastcall = .false.
256! **** needed but hardcoded
257
258elaps = (float(itimestep)-1.)*dt  ! elapsed seconds of simulation
259!----------------------------------------------!
260! what is done at the first step of simulation !
261!----------------------------------------------!
262IF (elaps .eq. 0.) THEN
263   flag_first_restart = .false.
264ELSE
265   flag_first_restart=flag_first_restart.OR.(.NOT. ALLOCATED(dp_save))
266ENDIF
267
268is_first_step: IF ((elaps .eq. 0.).or.flag_first_restart) THEN
269firstcall=.true.  !! for continuity with GCM, physics are always called at the first WRF timestep
270  !firstcall=.false. !! just in case you'd want to get rid of the physics
271test=0
272IF( .NOT. ALLOCATED( dp_save  ) ) THEN
273ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
274ALLOCATE(du_save(ngrid,nlayer))
275ALLOCATE(dv_save(ngrid,nlayer))
276ALLOCATE(dt_save(ngrid,nlayer))
277ALLOCATE(dq_save(ngrid,nlayer,nq))
278ENDIF
279dp_save(:)=0.    !! initialize these arrays ...
280du_save(:,:)=0.
281dv_save(:,:)=0.
282dt_save(:,:)=0.
283dq_save(:,:,:)=0.
284flag_first_restart=.false.
285
286ELSE ! is_first_step
287!-------------------------------------------------!
288! what is done for the other steps of simulation  !
289!-------------------------------------------------!
290IF (wappel_phys .gt. 0.) THEN
291firstcall = .false.
292test = MODULO(itimestep-1,int(wappel_phys)) ! WRF time is integrated time (itimestep=1 at the end of first step)
293                                            ! -- same strategy as diagfi in the LMD GCM
294                                            ! -- arguments of modulo must be of the same kind (here: integers)
295ELSE
296firstcall = .false.
297test = 9999   
298ENDIF
299ENDIF is_first_step
300
301!------------------------------------!
302! print info about domain initially  !
303! ... and whenever domain is changed !
304!------------------------------------!
305IF (previous_id .ne. id) THEN
306    print *,'** ',planet_type,' ** DOMAIN',id
307    print *, '** ',planet_type,' ** ... INITIALIZE DOMAIN',id
308    print *, '** ',planet_type,' ** ... PREVIOUS DOMAIN was',previous_id
309    print *, 'ITIMESTEP: ',itimestep
310    print *, 'TILES: ', i_start,i_end, j_start, j_end  ! numbers for simple runs, arrays for parallel runs
311    print *, 'DOMAIN: ', ids, ide, jds, jde, kds, kde
312    print *, 'MEMORY: ', ims, ime, jms, jme, kms, kme
313    print *, 'COMPUT: ', ips, ipe, jps, jpe, kps, kpe
314    print *, 'SIZE OF PHYSICS GRID for this process: ',ngrid
315    print *, 'ADVECTED TRACERS: ', num_3d_s-1
316    print *, 'PHYSICS IS CALLED EACH...',wappel_phys
317    print *, '-- means: PHYSICAL TIMESTEP=',zdt_split
318ENDIF
319
320
321!!!***********!!
322!!! IDEALIZED !!
323!!!***********!!
324IF (.not.(JULYR .le. 8999)) THEN
325  IF (firstcall .EQV. .true.) THEN
326    PRINT *,'** ',planet_type,'** IDEALIZED SIMULATION'
327    PRINT *,'** ',planet_type,'** BEWARE: input_coord must be here'
328  ENDIF
329  open(unit=14,file='input_coord',form='formatted',status='old')
330  rewind(14)
331  read(14,*) lon_input
332  read(14,*) lat_input
333  read(14,*) ls_input
334  read(14,*) lct_input
335  close(14)
336ENDIF
337
338!----------!
339! ALLOCATE !
340!----------!
341CALL allocate_interface(ngrid,nlayer,nq)
342!!!
343!!! BIG LOOP : 1. no call for physics, used saved values
344!!!
345IF (test.NE.0) THEN
346print *,'** ',planet_type,'** NO CALL FOR PHYSICS, go to next step...',test
347print*,'else'
348zdpsrf_omp(:)=dp_save(:)
349zdufi_omp(:,:)=du_save(:,:)
350zdvfi_omp(:,:)=dv_save(:,:)
351zdtfi_omp(:,:)=dt_save(:,:)
352zdqfi_omp(:,:,:)=dq_save(:,:,:)
353!!!
354!!! BIG LOOP : 2. calculate physical tendencies
355!!!
356ELSE ! if (test.EQ.0)
357!----------!
358! ALLOCATE !
359!----------!
360! interm
361ALLOCATE(dz8w_prof(nlayer))
362ALLOCATE(p8w_prof(nlayer+1))
363ALLOCATE(p_prof(nlayer))
364ALLOCATE(t_prof(nlayer))
365ALLOCATE(t8w_prof(nlayer))
366ALLOCATE(u_prof(nlayer))
367ALLOCATE(v_prof(nlayer))
368ALLOCATE(z_prof(nlayer))
369ALLOCATE(q_prof(nlayer,nq))
370
371
372!!!!!!!!!!!!!!!!!!!!!!!!!!!!
373!!!!!!!!!!!!!!!!!!!!!!!!!!!!
374!! PREPARE PHYSICS INPUTS !! 
375!!!!!!!!!!!!!!!!!!!!!!!!!!!!
376!!!!!!!!!!!!!!!!!!!!!!!!!!!!
377
378!! INITIALIZE AND ALLOCATE EVERYTHING !! here, only firstcall
379allocation_firstcall: IF (firstcall .EQV. .true.) THEN
380  !! PHYSICS VARIABLES (cf. iniphysiq in LMD GCM)
381  !! parameters are defined in the module_model_constants.F WRF routine
382  PRINT *,'** ',planet_type,'** INITIALIZE ARRAYS FOR PHYSICS'
383  !! need to get initial time first
384  CALL update_inputs_physiq_time(&
385            JULYR,JULDAY,GMT,&
386            elaps,&
387            lct_input,lon_input,ls_input,&
388            MY)
389  !! Set up initial time
390  phour_ini = JH_cur_split
391  !! Fill planetary parameters in modules
392  !! Values defined in the module_model_constants.F WRF routine
393  CALL update_inputs_physiq_constants
394  !! JL21 it seems nothing is done in update_inputs_physiq_constants for generic case !!!
395  !! Initialize physics
396  CALL iniphysiq(ngrid,nlayer,nq,wappel_phys,&
397                     wdaysec,floor(JD_cur), &
398                     1./reradius,g,r_d,cp,1)
399                     
400                     
401ENDIF allocation_firstcall
402
403!!*****************************!!
404!! PREPARE "FOOD" FOR PHYSIQ.F !!
405!!*****************************!!
406
407DO j = jps,jpe
408DO i = ips,ipe
409
410!!*******************************************!!
411!! FROM 3D WRF FIELDS TO 1D PHYSICS PROFILES !!
412!!*******************************************!!
413
414!-----------------------------------!
415! 1D subscript for physics "cursor" !
416!-----------------------------------!
417subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
418
419!--------------------------------------!
420! 1D columns : inputs for the physics  !
421! ... vert. coord., meteor. fields ... !
422!--------------------------------------!
423dz8w_prof(:) = dz8w(i,kps:kpe,j)   ! dz between full levels (m)   
424p_prof(:) = p_hyd(i,kps:kpe,j)         ! hydrostatic pressure at layers >> zplay_omp
425p8w_prof(:) = p_hyd_w(i,kps:kpe+1,j)         ! hydrostatic pressure at levels
426!JL22 using hydrostatic pressures to conserve mass
427t_prof(:) = t(i,kps:kpe,j)         ! temperature half level (K) >> pt
428t8w_prof(:) = t8w(i,kps:kpe,j)     ! temperature full level (K)
429u_prof(:) = u(i,kps:kpe,j)         ! zonal wind (A-grid: unstaggered) half level (m/s) >> zufi_omp 
430v_prof(:) = v(i,kps:kpe,j)         ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv
431z_prof(:) = z(i,kps:kpe,j)         ! geopotential height half level (m) >> zphi_omp/g
432
433
434!--------------------------------!
435! specific treatment for tracers !
436!--------------------------------!
437IF (TRACER_MODE == 1) THEN
438    ! to be clean we should have an automatized process that makes sure that moist is sent to igcm_h2o_vap and etc.
439    q_prof(:,1) = SCALAR(i,kps:kpe,j,P_QH2O) / (1.d0 + SCALAR(i,kps:kpe,j,P_QH2O)) !! P_xxx is the index for variable xxx.
440    q_prof(:,2) = SCALAR(i,kps:kpe,j,P_QH2O_ICE) / (1.d0 + SCALAR(i,kps:kpe,j,P_QH2O))
441    ! conversion from mass mixing ratio in WRF to specific concentration in Physiq
442ELSE IF ((TRACER_MODE >= 42).AND.(TRACER_MODE <= 45)) THEN
443    ! to be clean we should have an automatized process that makes sure that moist is sent to igcm_h2o_vap and etc.
444    q_prof(:,1) = moist(i,kps:kpe,j,P_QV) / (1.d0 + moist(i,kps:kpe,j,P_QV)) !! P_xxx is the index for variable xxx.
445    q_prof(:,2) = moist(i,kps:kpe,j,P_QC) / (1.d0 + moist(i,kps:kpe,j,P_QV))
446    ! conversion from mass mixing ratio in WRF to specific concentration in Physiq
447ELSE
448    q_prof(:,1:nq) = SCALAR(i,kps:kpe,j,2:nq+1)  !! the names were set above !! one dummy tracer in WRF
449    !JL22 cannot normalize to moist here as we do not know if it has been initialized.
450ENDIF
451
452IF (firstcall .EQV. .true.) THEN
453  !-----------------!
454  ! Optional output !
455  !-----------------!
456  IF ( (i == ips) .AND. (j == jps) ) THEN
457  PRINT *,'z_prof ',z_prof
458  PRINT *,'dz8w_prof',dz8w_prof
459  PRINT *,'p8w_prof ',p8w_prof
460  PRINT *,'p_prof ',p_prof
461  PRINT *,'t_prof ',t_prof
462  PRINT *,'t8w_prof ',t8w_prof
463  PRINT *,'u_prof ',u_prof
464  PRINT *,'v_prof ',v_prof
465  PRINT *,'q_prof ',q_prof
466  ENDIF
467ENDIF
468
469!---------------------------------------------!
470! in LMD physics, geopotential must be        !
471! expressed with respect to the local surface !
472!---------------------------------------------!
473zphi_omp(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) )
474
475!--------------------------------!
476! Dynamic fields for LMD physics !
477!--------------------------------!
478zplev_omp(subs,1:nlayer+1) = p8w_prof(1:nlayer+1)  !! NB: last level: no data
479zplay_omp(subs,:) = p_prof(:)
480ztfi_omp(subs,:) = t_prof(:)
481zufi_omp(subs,:) = u_prof(:)
482zvfi_omp(subs,:) = v_prof(:)
483flxwfi_omp(subs,:) = 0   !! NB: not used in the physics, only diagnostic...
484!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485!! for IDEALIZED CASES ONLY
486!IF (.not.(JULYR .le. 8999)) zplev_omp(subs,nlayer+1)=0.  !! zplev_omp(subs,nlayer+1)=ptop >> NO !
487!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
488
489! NOTE:
490! IF ( zplev_omp(subs,nlayer+1) .le. 0 ) zplev_omp(subs,nlayer+1)=ptop
491! cree des diagnostics delirants et aleatoires dans le transfert radiatif
492
493!---------!
494! Tracers ! 
495!---------!
496zqfi_omp(subs,:,:) = q_prof(:,:)  !! traceurs generiques, seuls noms sont specifiques
497
498ENDDO
499ENDDO
500
501!---------------------------------------------------------!
502! Ground geopotential                                     !
503! (=g*HT(i,j), only used in the microphysics: newcondens) !
504!---------------------------------------------------------!
505phisfi_val=g*(z_prof(1)-dz8w_prof(1)/2.)   !! NB: z_prof are half levels, so z_prof(1) is not the surface
506
507!!*****************!!
508!! CLEAN THE PLACE !!
509!!*****************!!
510DEALLOCATE(dz8w_prof)
511DEALLOCATE(z_prof)
512DEALLOCATE(p8w_prof)
513DEALLOCATE(p_prof)
514DEALLOCATE(t_prof)
515DEALLOCATE(t8w_prof)
516DEALLOCATE(u_prof)
517DEALLOCATE(v_prof)
518DEALLOCATE(q_prof)
519
520
521!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
522!!! ONE DOMAIN CASE
523!!! --> we simply need to initialize static and physics inputs
524!!! SEVERAL DOMAINS
525!!! --> we update static and physics inputs each time
526!!!     to account for domain change
527!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
528pass_interface: IF ( (firstcall .EQV. .true.) .or. (max_dom .GT. 1) ) THEN
529PRINT *,'** ',planet_type,'** PASS INTERFACE. EITHER FIRSTCALL or NESTED SIMULATION.'
530!!*******************************************!!
531!!*******************************************!!
532CALL update_inputs_physiq_geom( &
533            ims,ime,jms,jme,&
534            ips,ipe,jps,jpe,&
535            JULYR,ngrid,nlayer,&
536            DX,DY,MSFT,&
537            lat_input, lon_input,&
538            XLAT,XLONG)
539!!!
540CALL update_inputs_physiq_slope( &
541            ims,ime,jms,jme,&
542            ips,ipe,jps,jpe,&
543            JULYR,&
544            SLPX,SLPY)
545!!!
546print*, 'num_soil_layers, nsoil', num_soil_layers, nsoil
547CALL update_inputs_physiq_soil( &
548            ims,ime,jms,jme,&
549            ips,ipe,jps,jpe,&
550            JULYR,nsoil,&
551            P_TI,CST_TI,&
552            P_ISOIL,P_DSOIL,&
553            P_TSOIL,P_TSURF)
554!!!
555CALL update_inputs_physiq_surf( &
556            ims,ime,jms,jme,&
557            ips,ipe,jps,jpe,&
558            JULYR,TRACER_MODE,&
559            P_ALBEDO,CST_AL,&
560            P_TSURF,P_EMISS,P_CO2ICE,&
561            P_GW,P_Z0,CST_Z0,&
562            P_H2OICE,&
563            phisfi_val)
564!!!
565CALL update_inputs_physiq_turb( &
566            ims,ime,jms,jme,kms,kme,&
567            ips,ipe,jps,jpe,&
568            RESTART,isles,&
569            P_Q2,P_WSTAR)
570!!!
571CALL update_inputs_physiq_rad( &
572            ims,ime,jms,jme,&
573            ips,ipe,jps,jpe,&
574            RESTART,&
575            P_FLUXRAD)
576!!!
577ENDIF pass_interface
578!!*******************************************!!
579!!*******************************************!!
580
581!!!!!!!!!!!!!!!!!!!!!!
582!!!!!!!!!!!!!!!!!!!!!!
583!! CALL LMD PHYSICS !! 
584!!!!!!!!!!!!!!!!!!!!!!
585!!!!!!!!!!!!!!!!!!!!!!
586
587!!! initialize tendencies (security)
588zdpsrf_omp(:)=0.
589zdufi_omp(:,:)=0.
590zdvfi_omp(:,:)=0.
591zdtfi_omp(:,:)=0.
592zdqfi_omp(:,:,:)=0.
593
594call_physics : IF (wappel_phys .ne. 0.) THEN
595!print *, '** ',planet_type,'** CALL TO LMD PHYSICS'
596!!!
597
598CALL update_inputs_physiq_time(&
599            JULYR,JULDAY,GMT,&
600            elaps,&
601            lct_input,lon_input,ls_input,&
602            MY)
603!!!
604CALL call_physiq(planet_type,ngrid,nlayer,nq,       &
605                       firstcall,lastcall)
606!!!
607
608!---------------------------------------------------------------------------------!
609! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
610!---------------------------------------------------------------------------------!
611dp_save(:)=zdpsrf_omp(:)
612du_save(:,:)=zdufi_omp(:,:)
613dv_save(:,:)=zdvfi_omp(:,:)
614dt_save(:,:)=zdtfi_omp(:,:)
615dq_save(:,:,:)=zdqfi_omp(:,:,:)
616
617!! OUTPUT OUTPUT OUTPUT
618!-------------------------------------------------------!
619! Save key variables for restart and output and nesting ! 
620!-------------------------------------------------------!
621!!!
622CALL update_outputs_physiq_surf( &
623            ims,ime,jms,jme,&
624            ips,ipe,jps,jpe,&
625            TRACER_MODE,&
626            P_TSURF,P_CO2ICE,&
627            P_H2OICE)
628!!!
629CALL update_outputs_physiq_soil( &
630            ims,ime,jms,jme,&
631            ips,ipe,jps,jpe,&
632            nsoil,&
633            P_TSOIL)
634!!!
635CALL update_outputs_physiq_rad( &
636            ims,ime,jms,jme,&
637            ips,ipe,jps,jpe,&
638            P_FLUXRAD)
639!!!
640CALL update_outputs_physiq_turb( &
641            ims,ime,jms,jme,kms,kme,&
642            ips,ipe,jps,jpe,kps,kpe,&
643            P_Q2,P_WSTAR,&
644            HFMAX,ZMAX,USTM,HFX)
645!!!
646CALL update_outputs_physiq_diag( &
647            ims,ime,jms,jme,kms,kme,&
648            ips,ipe,jps,jpe,kps,kpe,&
649            HR_SW,HR_LW,HR_DYN,DT_RAD,&
650            CLOUDFRAC,TOTCLOUDFRAC,RH,&
651            DQICE,DQVAP,REEVAP,SURFRAIN,&
652            ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
653            FLUXSURF_LW,FLXGRD,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF,LATENT_HF)
654!!!
655!print *, '** ',planet_type,'** OUTPUT PHYSICS DONE'
656
657ENDIF call_physics
658
659ENDIF    ! end of BIG LOOP 2.
660!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
661!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
662!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
663!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
664
665!!***************************!!
666!! DEDUCE TENDENCIES FOR WRF !!
667!!***************************!!
668RTHPLATEN(ims:ime,kms:kme,jms:jme)=0.
669RUPLATEN(ims:ime,kms:kme,jms:jme)=0.
670RVPLATEN(ims:ime,kms:kme,jms:jme)=0.
671PSFC(ims:ime,jms:jme)=p8w(ims:ime,kms,jms:jme) ! was done in surface driver in regular WRF
672!------------------------------------------------------------------!
673! WRF adds the physical and dynamical tendencies                   !
674! thus the tendencies must be passed as 'per second' tendencies    !
675! --when physics is not called, the applied physical tendency      !
676! --is the one calculated during the last call to physics          !
677!------------------------------------------------------------------!
678!print*,'pdt',pdt(1,1),pdt(1,nlayer)
679!print*,'exner',exner(1,:,1)
680DO j = jps,jpe
681DO i = ips,ipe
682
683  subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
684
685    ! zonal wind
686  RUPLATEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe)
687    ! meridional wind
688  RVPLATEN(i,kps:kpe,j) = zdvfi_omp(subs,kps:kpe)
689    ! potential temperature
690    ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
691  RTHPLATEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j)
692    ! update surface pressure (cf CO2 cycle in physics)
693    ! here dt is needed
694  PSFC(i,j)=PSFC(i,j)+zdpsrf_omp(subs)*dt
695    ! tracers
696  SCALAR(i,kps:kpe,j,1)=0.
697  SELECT CASE (TRACER_MODE)
698    CASE(0)
699      SCALAR(i,kps:kpe,j,:)=0.
700    CASE(1)
701      scalar(i,kps:kpe,j,P_QH2O)=scalar(i,kps:kpe,j,P_QH2O) &
702          +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+scalar(i,kps:kpe,j,P_QH2O))
703      scalar(i,kps:kpe,j,P_QH2O_ICE)=scalar(i,kps:kpe,j,P_QH2O_ICE) &
704           +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+scalar(i,kps:kpe,j,P_QH2O))
705       ! if you want to use this mode, RTHPLATEN should be corrected as below.
706       ! we keep it like that for the moment for testing.
707    CASE(42)
708      moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) &
709          +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
710      moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) &
711           +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
712      RTHPLATEN(i,kps:kpe,j) = RTHPLATEN(i,kps:kpe,j) &
713           * (1.d0+moist(i,kps:kpe,j,P_QV))/(1.d0+rvovrd*moist(i,kps:kpe,j,P_QV))
714           ! correct dT/dt assuming a constant molar heat capacity.
715           ! Specific heat cappacity scales with molar mass.
716    CASE(43)
717      moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) &
718          +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
719      moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) &
720           +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
721      tau_decay=86400.*100. !! why not make it a namelist argument?
722      SCALAR(i,kps:kpe,j,P_MARKER) = SCALAR(i,kps:kpe,j,P_MARKER)*exp(-dt/tau_decay)
723      SCALAR(i,1,j,P_MARKER) = 1. !! this tracer is emitted in the surface layer
724      RTHPLATEN(i,kps:kpe,j) = RTHPLATEN(i,kps:kpe,j) &
725        * (1.d0+moist(i,kps:kpe,j,P_QV))/(1.d0+rvovrd*moist(i,kps:kpe,j,P_QV))
726        ! correct dT/dt assuming a constant molar heat capacity.
727        ! Specific heat cappacity scales with molar mass.
728    CASE(44)
729      moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) &
730          +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
731      moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) &
732           +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
733    CASE(45)
734      moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) &
735          +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
736      moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) &
737           +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
738      tau_decay=86400.*100. !! why not make it a namelist argument?
739      SCALAR(i,kps:kpe,j,P_MARKER) = SCALAR(i,kps:kpe,j,P_MARKER)*exp(-dt/tau_decay)
740      SCALAR(i,1,j,P_MARKER) = 1. !! this tracer is emitted in the surface layer
741    CASE DEFAULT
742      !SCALAR(i,kps:kpe,j,2:nq+1)=SCALAR(i,kps:kpe,j,2:nq+1)+zdqfi_omp(subs,kps:kpe,1:nq)*dt !!! here dt is needed
743      scalar(i,kps:kpe,j,2:nq+1)=scalar(i,kps:kpe,j,2:nq+1) &
744           +zdqfi_omp(subs,kps:kpe,1:nq)*dt
745  END SELECT
746
747ENDDO
748ENDDO
749CALL deallocate_interface
750!!*****!!
751!! END !!
752!!*****!!
753!print *,'** ',planet_type,'** END LMD PHYSICS'
754previous_id = id
755!----------------------------------------------------------------!
756! use debug (see solve_em) if you wanna display some message ... !
757!----------------------------------------------------------------!
758END SUBROUTINE lmd_driver
759
760END MODULE module_lmd_driver
Note: See TracBrowser for help on using the repository browser.