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

Last change on this file since 3661 was 3661, checked in by emoisan, 4 months ago

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