source: trunk/WRF.COMMON/INTERFACES/module_lmd_driver.F @ 3529

Last change on this file since 3529 was 2758, checked in by aslmd, 2 years ago

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