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

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

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