source: trunk/WRF.COMMON/WRFV2/phys/module_lmd_driver.F.old @ 3567

Last change on this file since 3567 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: 54.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!*******************************************************************************
15MODULE module_lmd_driver
16CONTAINS
17    SUBROUTINE lmd_driver(id,max_dom,DT,ITIMESTEP,XLAT,XLONG, &
18        IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME,KMS,KME, &
19        i_start,i_end,j_start,j_end,kts,kte,num_tiles, &
20        DX,DY, &
21        MSFT,MSFU,MSFV, &
22        GMT,JULYR,JULDAY, &
23        P8W,DZ8W,T8W,Z,HT, &
24        U,V,W,TH,T,P,EXNER,RHO, &
25        PTOP, &
26        RADT,CUDT, &
27        TSK,PSFC, &
28        RTHBLTEN,RUBLTEN,RVBLTEN, &
29        num_3d_s,SCALAR, &
30        num_3d_m,MOIST,&
31        MARS_MODE, &
32        M_ALBEDO,M_TI,M_CO2ICE,M_EMISS, &
33        M_H2OICE, &
34        M_TSOIL, &
35        M_Q2, &
36        M_TSURF, &
37#ifdef NEWPHYS
38        M_FLUXRAD, &
39        M_WSTAR, &
40        M_ISOIL, &
41        M_DSOIL, &
42        M_Z0, &
43        CST_Z0, &
44#endif
45        M_GW, &
46        NUM_SOIL_LAYERS, &
47        CST_AL, &
48        CST_TI, &
49        isfflx, &
50        diff_opt, &
51        km_opt, &
52        HISTORY_INTERVAL, &
53#ifndef NOPHYS
54#include "module_lmd_driver_output1.inc"
55#endif
56        SLPX,SLPY,RESTART)
57! NB: module_lmd_driver_output1.inc : output arguments generated from Registry
58
59
60
61
62
63!==================================================================
64! USES
65!==================================================================
66   USE module_model_constants
67   USE module_wrf_error
68! add new modules here, if needed ...
69
70!==================================================================
71  IMPLICIT NONE
72!==================================================================
73
74!==================================================================
75! COMMON
76!==================================================================
77
78#ifndef NOPHYS
79!! the only common needed is the one defining the physical grid
80!! -- path is hardcoded, but the structure is not subject to change
81!! -- please put # if needed by the pre-compilation process     
82!
83!
84include "../mars_lmd/libf/phymars/dimphys.h"
85!
86!--to be commented because there are tests in the physics ?
87!--TODO: get rid of the ...mx first in this routine and .inc
88
89!
90! INCLUDE AUTOMATIQUEMENT GENERE A PARTIR DU REGISTRY
91!
92include "../mars_lmd/libf/phymars/wrf_output_2d.h"
93include "../mars_lmd/libf/phymars/wrf_output_3d.h"
94#endif
95
96!==================================================================
97! VARIABLES
98!==================================================================
99
100! WRF Dimensions
101INTEGER, INTENT(IN   )    ::   &
102      ids,ide,jds,jde,kds,kde, &
103      ims,ime,jms,jme,kms,kme, &
104            kts,kte,num_tiles, &
105              NUM_SOIL_LAYERS
106INTEGER, DIMENSION(num_tiles), INTENT(IN) ::  &
107      i_start,i_end,j_start,j_end   
108! Scalars
109INTEGER, INTENT(IN  ) :: JULDAY, itimestep, julyr,id,max_dom,MARS_MODE
110INTEGER, INTENT(IN  ) :: isfflx,diff_opt,km_opt
111REAL, INTENT(IN  ) :: GMT,dt,dx,dy,RADT,CUDT
112REAL, INTENT(IN  ) :: CST_AL, CST_TI
113REAL, INTENT(IN  ) :: PTOP
114! 2D arrays         
115REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  :: &
116     MSFT,MSFU,MSFV, &
117     XLAT,XLONG,HT,  &
118     M_ALBEDO,M_TI,M_EMISS, &
119     SLPX,SLPY
120REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT   )  :: &
121     M_CO2ICE,M_H2OICE, &
122     M_TSURF
123! 3D arrays
124REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
125     dz8w,p8w,p,exner,t,t8w,rho,u,v,w,z,th
126!REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: &
127!     M_Q2
128REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: &
129     M_Q2
130!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
131INTEGER, INTENT(IN   ) :: HISTORY_INTERVAL
132!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
133REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(INOUT )  :: &
134     M_TSOIL
135#ifdef NEWPHYS
136REAL, INTENT(IN  ) :: CST_Z0
137REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(IN   )  :: &
138     M_ISOIL, M_DSOIL         
139REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  :: &
140     M_Z0
141REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT   )  :: &
142     M_FLUXRAD,M_WSTAR
143#endif
144REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: &
145     M_GW
146! 4D arrays
147INTEGER, INTENT(IN ) :: num_3d_s,num_3d_m
148REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_s ), INTENT(INOUT ) :: &
149     scalar
150REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_m ), INTENT(INOUT ) :: &
151     moist
152! Logical
153LOGICAL, INTENT(IN ) :: restart
154
155!-------------------------------------------
156! OUTPUT VARIABLES
157!-------------------------------------------
158!
159! Generated from Registry
160!
161! default definitions :
162! 2D : TSK, PSFC
163! 3D : RTHBLTEN,RUBLTEN,RVBLTEN
164#ifndef NOPHYS
165#include "module_lmd_driver_output2.inc"
166   REAL, DIMENSION(:,:), ALLOCATABLE :: output_tab2d
167   REAL, DIMENSION(:,:,:), ALLOCATABLE :: output_tab3d
168#else
169   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: PSFC,TSK
170   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT)  :: RTHBLTEN,RUBLTEN,RVBLTEN
171   REAL,DIMENSION(:),ALLOCATABLE,SAVE :: dtrad
172#endif
173!-------------------------------------------
174! OUTPUT VARIABLES
175!-------------------------------------------
176
177
178!-------------------------------------------   
179! LOCAL  VARIABLES
180!-------------------------------------------
181   INTEGER ::    i,j,k,its,ite,jts,jte,ij,kk
182   INTEGER ::    subs,iii
183
184   ! *** for Mars Mode 20 and 21 ***
185   REAL ::    tau_decay, lct, zilctmin
186   INTEGER, SAVE :: zipbl(500) !index of zi in file input_zipbl
187   REAL, SAVE :: zilct(500) !corresponding local time in input_zipbl
188   INTEGER :: lctindex,ziindex
189   LOGICAL :: end_of_file
190   ! *** ----------------------- ***
191
192   ! *** for LMD physics
193   ! ------> inputs:
194   INTEGER :: ngrid,nlayer,nq,nsoil
195   REAL :: pday,ptime,MY 
196   REAL :: aire_val,lat_val,lon_val
197   REAL :: phisfi_val,albedodat_val,inertiedat_val
198   REAL :: tsurf_val,co2ice_val,emis_val
199   REAL :: zmea_val,zstd_val,zsig_val,zgam_val,zthe_val
200   REAL :: theta_val, psi_val
201   LOGICAL :: firstcall,lastcall,tracerdyn
202   REAL,DIMENSION(:),ALLOCATABLE :: q2_val, qsurf_val, tsoil_val
203#ifdef NEWPHYS
204   REAL :: wstar_val,fluxrad_val
205   REAL,DIMENSION(:),ALLOCATABLE :: isoil_val, dsoil_val
206   REAL :: z0_val
207#endif
208   REAL,DIMENSION(:),ALLOCATABLE :: aire_vec,lat_vec,lon_vec
209   REAL,DIMENSION(:),ALLOCATABLE :: walbedodat,winertiedat,wphisfi
210   REAL,DIMENSION(:),ALLOCATABLE :: wzmea,wzstd,wzsig,wzgam,wzthe
211   REAL,DIMENSION(:),ALLOCATABLE :: wtheta, wpsi
212   ! v--- can they be modified ?
213   REAL,DIMENSION(:),ALLOCATABLE :: wtsurf,wco2ice,wemis
214   REAL,DIMENSION(:,:),ALLOCATABLE :: wq2,wqsurf,wtsoil
215#ifdef NEWPHYS
216   REAL,DIMENSION(:),ALLOCATABLE :: wwstar,wfluxrad
217   REAL,DIMENSION(:),ALLOCATABLE :: wz0tab
218   REAL,DIMENSION(:,:),ALLOCATABLE :: wisoil,wdsoil
219   CHARACTER*20,DIMENSION(:),ALLOCATABLE :: wtnom
220#endif
221   ! ----------
222   REAL,DIMENSION(:,:),ALLOCATABLE :: pplev,pplay,pphi,pu,pv,pt,pw
223   REAL,DIMENSION(:,:,:),ALLOCATABLE :: pq 
224
225   ! <------ outputs:
226   !     physical tendencies
227   REAL,DIMENSION(:),ALLOCATABLE :: pdpsrf
228   REAL,DIMENSION(:,:),ALLOCATABLE :: pdu,pdv,pdt
229   REAL,DIMENSION(:,:,:),ALLOCATABLE :: pdq
230   ! ... intermediate arrays
231   REAL, DIMENSION(:), ALLOCATABLE  :: &
232     dz8w_prof,p8w_prof,p_prof,t_prof,t8w_prof, &
233     u_prof,v_prof,z_prof
234!!   pi_prof, rho_prof, th_prof, &
235#ifdef NEWPHYS
236   REAL, DIMENSION(:,:), ALLOCATABLE :: q_prof
237#else
238   REAL, DIMENSION(:), ALLOCATABLE  :: &
239     water_vapor_prof, water_ice_prof
240#endif
241
242
243   ! Additional control variables
244   INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe
245   REAL :: elaps, ptimestep
246   INTEGER :: wday_ini, test, test2
247   REAL :: wappel_phys
248   LOGICAL :: flag_LES
249   LOGICAL, SAVE :: flag_first_restart
250!**************************************************
251! IMPORTANT: pour travailler avec grid nesting
252! --- deux comportements distincts du save
253! ... ex1: ferme planeto avec PGF+MPI: standard
254! ... ex2: iDataPlex avec IFORT+MPI: SPECIAL_NEST_SAVE
255!**************************************************
256#ifdef SPECIAL_NEST_SAVE
257      !  une dimension supplementaire liee au nest
258      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
259             dp_save
260      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
261             du_save, dv_save, dt_save
262      REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
263             dq_save     
264#ifndef NORESTART
265      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
266             save_tsoil_restart
267      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
268             save_tsurf_restart
269      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
270             save_co2ice_restart
271      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
272             save_q2_restart
273      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
274             save_qsurf_restart
275#ifdef NEWPHYS
276      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
277             save_wstar_restart
278      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
279             save_fluxrad_restart
280#endif
281#endif
282#else
283      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
284             dp_save
285      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
286             du_save, dv_save, dt_save
287      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
288             dq_save     
289#ifndef NORESTART
290!! FOR RESTART
291      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
292             save_tsoil_restart
293      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
294             save_tsurf_restart
295      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
296             save_co2ice_restart
297      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
298             save_q2_restart
299      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
300             save_qsurf_restart
301#ifdef NEWPHYS
302      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
303             save_wstar_restart
304      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
305             save_fluxrad_restart
306#endif
307#endif
308#endif
309
310!!!IDEALIZED IDEALIZED
311      REAL :: lat_input, lon_input, ls_input, lct_input
312!!!IDEALIZED IDEALIZED
313
314!!!!!!!!!!!!!! TEST NaN or Inf : ne fonction qu'avec r4 i4
315integer :: mask_nan = 2139095040
316integer :: chuck_norris
317!!!!!!!!!!!!!! TEST NaN or Inf : ne fonction qu'avec r4 i4
318
319!REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
320!     UMAX, UMIN, VMAX, VMIN, WMAX, WMIN, TMAX, TMIN
321
322
323
324!==================================================================
325! CODE
326!==================================================================
327
328IF (JULYR .ne. 9999) THEN
329    flag_LES = .false.  ! "True" LES is not available in this version
330    PRINT *, '*** REAL-CASE SIMULATION ***'
331ELSE
332     PRINT *, '*** IDEALIZED SIMULATION ***'
333     IF ((diff_opt .eq. 2) .and. (km_opt .eq. 2)) THEN
334          PRINT *, '*** type: LES ***'
335          PRINT *, '*** diff_opt = 2 *** km_opt = 2'
336          PRINT *, '*** forcing is isfflx = ',isfflx
337          flag_LES = .true.
338          !! SPECIAL LES
339     ELSE
340          PRINT *, '*** type: not LES ***'
341          PRINT *, '*** diff_opt = ',diff_opt
342          PRINT *, '*** km_opt = ',km_opt
343          flag_LES = .false.
344          !! IDEALIZED, no LES
345          !! cependant, ne veut-on pas pouvoir
346          !! prescrire un flux en idealise ??
347     ENDIF
348ENDIF
349
350
351print *,'** Mars ** DOMAIN',id
352
353!-------------------------!
354! TWEAK on WRF DIMENSIONS !
355!-------------------------!
356its = i_start(1)   ! define here one big tile
357ite = i_end(num_tiles)
358jts = j_start(1)
359jte = j_end(num_tiles)
360!!
361IF (flag_LES .eqv. .false.) THEN
362 relax=0
363 sponge_top=0               ! another value than 0 triggers instabilities 
364 IF (id .gt. 1) relax=2     ! fix to avoid noise in nesting simulations ; 1 >> too much noise ...
365ENDIF
366ips=its
367ipe=ite
368jps=jts
369jpe=jte
370IF (flag_LES .eqv. .false.) THEN
371 IF (ips .eq. ids)   ips=its+relax !! IF tests necesary for parallel runs
372 IF (ipe .eq. ide-1) ipe=ite-relax
373 IF (jps .eq. jds)   jps=jts+relax
374 IF (jpe .eq. jde-1) jpe=jte-relax
375ENDIF
376kps=kts         !! start at surface
377IF (flag_LES .eqv. .false.) THEN
378 kpe=kte-sponge_top
379ELSE
380 PRINT *, '*** IDEALIZED SIMULATION: LES *** kpe=kte'
381 kpe=kte !-sponge_top
382ENDIF
383
384!----------------------------!
385! DIMENSIONS FOR THE PHYSICS !
386!----------------------------!
387wday_ini = JULDAY - 1      !! GCM convention
388wappel_phys = RADT
389ptimestep = dt*wappel_phys            ! physical timestep (s)
390ngrid=(ipe-ips+1)*(jpe-jps+1)         ! size of the horizontal grid
391nlayer = kpe-kps+1                    ! number of vertical layers: nlayermx
392nsoil = NUM_SOIL_LAYERS               ! number of soil layers: nsoilmx
393if (num_3d_s > 1) then                ! number of advected fields
394        nq = num_3d_s-1               
395else
396        nq = 1
397endif
398! **** needed but hardcoded
399lastcall = .false.
400! **** needed but hardcoded
401
402          !PRINT *, ips, ipe, jps, jpe
403          !PRINT *, ngrid
404
405elaps = (float(itimestep)-1.)*dt  ! elapsed seconds of simulation
406!----------------------------------------------!
407! what is done at the first step of simulation !
408!----------------------------------------------!
409IF (elaps .eq. 0.) THEN
410   flag_first_restart = .false.
411ELSE
412   flag_first_restart=flag_first_restart.OR.(.NOT. ALLOCATED(dp_save))
413ENDIF
414
415IF ((elaps .eq. 0.).or.flag_first_restart) THEN
416firstcall=.true.  !! for continuity with GCM, physics are always called at the first WRF timestep
417  !firstcall=.false. !! just in case you'd want to get rid of the physics
418test=0
419#ifdef SPECIAL_NEST_SAVE
420PRINT *, 'ALLOCATED dp_save ???', ALLOCATED( dp_save ), id
421IF( .NOT. ALLOCATED( dp_save  ) ) THEN
422   PRINT *, '**** check **** OK I ALLOCATE one save ARRAY for all NESTS ', max_dom, id
423   !! here are the arrays that would be useful to save physics tendencies
424   ALLOCATE(dp_save(ngrid,max_dom))
425   ALLOCATE(du_save(ngrid,nlayer,max_dom))
426   ALLOCATE(dv_save(ngrid,nlayer,max_dom))
427   ALLOCATE(dt_save(ngrid,nlayer,max_dom))
428   ALLOCATE(dq_save(ngrid,nlayer,nq,max_dom))
429   dp_save(:,:)=0.    !! initialize these arrays ...
430   du_save(:,:,:)=0.
431   dv_save(:,:,:)=0.
432   dt_save(:,:,:)=0.
433   dq_save(:,:,:,:)=0.
434#ifndef NORESTART
435   ! Restart save arrays
436   ALLOCATE(save_tsoil_restart(ngrid,nsoil,max_dom))
437   ALLOCATE(save_co2ice_restart(ngrid,max_dom))
438   ALLOCATE(save_q2_restart(ngrid,nlayer+1,max_dom))
439   ALLOCATE(save_qsurf_restart(ngrid,nq,max_dom))
440   ALLOCATE(save_tsurf_restart(ngrid,max_dom))
441   save_tsoil_restart(:,:,:)=0.
442   save_co2ice_restart(:,:)=0.
443   save_q2_restart(:,:,:)=0.
444   save_qsurf_restart(:,:,:)=0.
445   save_tsurf_restart(:,:)=0.
446#ifdef NEWPHYS
447   ALLOCATE(save_wstar_restart(ngrid,max_dom))
448   ALLOCATE(save_fluxrad_restart(ngrid,max_dom))
449   save_wstar_restart(:,:)=0.
450   save_fluxrad_restart(:,:)=0.
451#endif
452#endif
453ENDIF
454IF (id .lt. max_dom) THEN
455   flag_first_restart=.true.
456ELSE
457   flag_first_restart=.false.
458ENDIF
459#else
460IF( .NOT. ALLOCATED( dp_save  ) ) THEN
461ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
462ALLOCATE(du_save(ngrid,nlayer))
463ALLOCATE(dv_save(ngrid,nlayer))
464ALLOCATE(dt_save(ngrid,nlayer))
465ALLOCATE(dq_save(ngrid,nlayer,nq))
466ENDIF
467dp_save(:)=0.    !! initialize these arrays ...
468du_save(:,:)=0.
469dv_save(:,:)=0.
470dt_save(:,:)=0.
471dq_save(:,:,:)=0.
472flag_first_restart=.false.
473#ifndef NORESTART
474! Restart save arrays
475IF( .NOT. ALLOCATED( save_tsoil_restart  ) ) THEN
476ALLOCATE(save_tsoil_restart(ngrid,nsoil))
477ALLOCATE(save_co2ice_restart(ngrid))
478ALLOCATE(save_q2_restart(ngrid,nlayer+1))
479ALLOCATE(save_qsurf_restart(ngrid,nq))
480ALLOCATE(save_tsurf_restart(ngrid))
481ENDIF
482save_tsoil_restart(:,:)=0.
483save_co2ice_restart(:)=0.
484save_q2_restart(:,:)=0.
485save_qsurf_restart(:,:)=0.
486save_tsurf_restart(:)=0.
487#ifdef NEWPHYS
488IF( .NOT. ALLOCATED( save_wstar_restart  ) ) THEN
489ALLOCATE(save_wstar_restart(ngrid))
490ALLOCATE(save_fluxrad_restart(ngrid))
491ENDIF
492save_wstar_restart(:)=0.
493save_fluxrad_restart(:)=0.
494#endif
495#endif
496#endif
497
498!! put here some general information you'd like to print just once
499    print *, 'TILES: ', i_start,i_end, j_start, j_end  ! numbers for simple runs, arrays for parallel runs
500    print *, 'DOMAIN: ', ids, ide, jds, jde
501    print *, 'MEMORY: ', ims, ime, jms, jme
502    print *, 'ADVECTED TRACERS: ', num_3d_s-1
503    print *, 'PHYSICS IS CALLED EACH...',wappel_phys
504!! put here some general information you'd like to print just once
505ELSE
506!-------------------------------------------------!
507! what is done for the other steps of simulation  !
508!-------------------------------------------------!
509IF (wappel_phys .gt. 0.) THEN
510firstcall = .false.
511test = MODULO(itimestep-1,int(wappel_phys)) ! WRF time is integrated time (itimestep=1 at the end of first step)
512                                            ! -- same strategy as diagfi in the LMD GCM
513                                            ! -- arguments of modulo must be of the same kind (here: integers)
514ELSE
515firstcall = .false.
516test = 9999   
517ENDIF
518ENDIF
519
520!!!!! for 'subgrid' temporal diagnostics
521!test2 = MODULO(elaps,history_interval*100.)
522
523!!!******!!
524!!! TIME !!
525!!!******!!
526IF (JULYR .ne. 9999) THEN
527  !
528  ! specified
529  !
530  ptime = (GMT + elaps/3700.) !! universal time (0<ptime<1): ptime=0.5 at 12:00 UT
531  ptime = MODULO(ptime,24.)   !! the two arguments of MODULO must be of the same type
532  ptime = ptime / 24.
533  pday = (JULDAY - 1 + INT((3700*GMT+elaps)/88800))
534  pday = MODULO(int(pday),669)
535  MY = (JULYR-2000) + (88800.*(JULDAY - 1)+3700.*GMT+elaps)/59496000.
536  MY = INT(MY)
537ELSE
538  !
539  ! idealized
540  !
541  PRINT *,'** Mars ** IDEALIZED SIMULATION'
542  PRINT *,'** Mars ** BEWARE: input_coord must be here'
543  open(unit=14,file='input_coord',form='formatted',status='old')
544  rewind(14)
545  read(14,*) lon_input
546  read(14,*) lat_input
547  read(14,*) ls_input
548  read(14,*) lct_input
549  close(14)
550  ptime = lct_input - lon_input / 15. + elaps/3700.
551  ptime = MODULO(ptime,24.)
552  ptime = ptime / 24.
553  pday = floor(ls2sol(ls_input)) + INT((3700*(lct_input - lon_input / 15.) + elaps)/88800)
554  pday = MODULO(int(pday),669)
555  MY = 2024
556  wday_ini = floor(ls2sol(ls_input))
557ENDIF
558print *,'** Mars ** TIME IS', pday, ptime*24.
559
560
561!!****************!!
562!! CHECK DYNAMICS !!
563!!****************!!
564!IF ((MAXVAL(t) > 500).OR.(MINVAL(t,MASK = t > 0) <= 50)) THEN
565!        PRINT *,'******************   CRASH   *******************'
566!        PRINT *,'Irrealistic temperature...', MAXLOC(t), MINLOC(t)
567!PRINT *, t
568!        PRINT *,'************************************************'
569!        STOP
570!ENDIF   
571
572!print *, 'check dynamics'
573!  !!! in some cases, weird values are displayed
574!  !!! despite the fact that outputs are OK...
575!  !print *, 'u', MAXVAL(u), MINVAL(u)
576!  !print *, 'v', MAXVAL(v), MINVAL(v)
577!  !print *, 'w', MAXVAL(w), MINVAL(w)
578!  !print *, 't', MAXVAL(t), MINVAL(t, MASK = t > 0)
579!print *, 'u', u(10,1,10), u(10,15,10)
580!print *, 'v', v(10,1,10), v(10,15,10)
581!print *, 'w', w(10,1,10), w(10,15,10)
582!print *, 't', t(10,1,10), t(10,15,10)
583
584!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
585!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
586!IF (test2.EQ.0) THEN
587!  print *, 'compute stats'
588!  print *, 'RESET'
589!   uave  = uave*0.
590!     vave  = vave*0.
591!     tave  = tave*0.
592!     wave  = wave*0.
593!   ustd  = ustd*0.
594!     vstd  = vstd*0.
595!     tstd  = tstd*0.
596!     wstd  = wstd*0.
597!ENDIF
598!  uave = uave + u  * dt / (float(history_interval)*100.)
599!  vave = vave + v  * dt / (float(history_interval)*100.)
600!  tave = tave + th * dt / (float(history_interval)*100.)
601!  wave = wave + w  * dt / (float(history_interval)*100.)
602!  ustd = ustd + u  * u  * dt / (float(history_interval)*100.)
603!  vstd = vstd + v  * v  * dt / (float(history_interval)*100.)
604!  tstd = tstd + th * th * dt / (float(history_interval)*100.)
605!  wstd = wstd + w  * w  * dt / (float(history_interval)*100.)
606!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
607!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
608
609!!!!!
610!!!!! PROBLEME avec 3 NESTS !!!
611!!!!!
612!chuck_norris = transfer( u(1,1,1),chuck_norris )        !! astuce J. Lefrere - test NaN or Inf
613!PRINT *,'check stability'
614!IF ( iand( chuck_norris,mask_nan ) == mask_nan ) THEN   !! astuce J. Lefrere - test NaN or Inf
615!        PRINT *, u(1,1,1)
616!        PRINT *,'******************   CRASH   *******************'
617!        PRINT *,'************************************************'
618!        PRINT *,'NaN or Inf appeared in the simulation ...'
619!        PRINT *,'...this may be due to numerical or dynamical instability'
620!        PRINT *,'************************************************'
621!        PRINT *,'POSSIBLE SOLUTIONS:'
622!        PRINT *,'>> IF nonhydrostatic mode,'
623!        PRINT *,'   --> check that smdiv, emdiv and epssm are not 0.'
624!        PRINT *,'>> IF cfl is violated, '
625!        PRINT *,'   --> try to lower the dynamical timestep'
626!        PRINT *,'>> IF topographical gradients are high near specified bdy,'
627!        PRINT *,'   --> try to redefine the domain'
628!        PRINT *,'************************************************'
629!        STOP
630!ELSE
631!        PRINT *,'OK OK OK OK'
632!ENDIF
633!ENDIF
634        !IF (          ANY(isNaN(u)) &
635        !         .OR. ANY(isNaN(v)) &
636        !         .OR. ANY(isNaN(t)) ) THEN
637        ! >>> ne marche qu'avec g95
638!print *, 'check dynamics'
639!print *, 'u', MAXVAL(u), MINVAL(u)
640!print *, 'v', MAXVAL(v), MINVAL(v)
641!print *, 't', MAXVAL(t), MINVAL(t, MASK = t > 0)
642
643
644
645!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
646!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
647!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
648!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
649!----------!
650! ALLOCATE !
651!----------!
652ALLOCATE(pdpsrf(ngrid))
653ALLOCATE(pdu(ngrid,nlayer))
654ALLOCATE(pdv(ngrid,nlayer))
655ALLOCATE(pdt(ngrid,nlayer))
656ALLOCATE(pdq(ngrid,nlayer,nq))
657!!!
658!!! BIG LOOP : 1. no call for physics, used saved values
659!!!
660IF (test.NE.0) THEN
661print *,'** Mars ** NO CALL FOR PHYSICS, go to next step...',test
662#ifdef SPECIAL_NEST_SAVE
663pdpsrf(:)=dp_save(:,id)
664pdu(:,:)=du_save(:,:,id)
665pdv(:,:)=dv_save(:,:,id)
666pdt(:,:)=dt_save(:,:,id)
667pdq(:,:,:)=dq_save(:,:,:,id)
668#else
669pdpsrf(:)=dp_save(:)
670pdu(:,:)=du_save(:,:)
671pdv(:,:)=dv_save(:,:)
672pdt(:,:)=dt_save(:,:)
673pdq(:,:,:)=dq_save(:,:,:)
674#endif
675!!!
676!!! BIG LOOP : 2. calculate physical tendencies
677!!!
678ELSE
679!----------!
680! ALLOCATE !
681!----------!
682! inputs ...
683IF (firstcall .EQV. .true.) THEN
684  ALLOCATE(aire_vec(ngrid))!
685  ALLOCATE(lon_vec(ngrid))!
686  ALLOCATE(lat_vec(ngrid))!
687  ALLOCATE(walbedodat(ngrid))!
688  ALLOCATE(winertiedat(ngrid))!
689  ALLOCATE(wphisfi(ngrid))!
690  ALLOCATE(wzmea(ngrid))!
691  ALLOCATE(wzstd(ngrid))!
692  ALLOCATE(wzsig(ngrid))!
693  ALLOCATE(wzgam(ngrid))!
694  ALLOCATE(wzthe(ngrid))!
695  ALLOCATE(wtheta(ngrid))!
696  ALLOCATE(wpsi(ngrid))!
697#ifdef NEWPHYS
698  ALLOCATE(wz0tab(ngrid))!
699#endif
700ENDIF
701ALLOCATE(wtsurf(ngrid))          !!!!!
702#ifndef NOPHYS
703ALLOCATE(output_tab2d(ngrid,n2d))
704ALLOCATE(output_tab3d(ngrid,nlayer,n3d))
705#endif
706ALLOCATE(wco2ice(ngrid))         !!!!!
707ALLOCATE(wemis(ngrid))           !!!!!
708ALLOCATE(q2_val(nlayer+1))
709ALLOCATE(qsurf_val(nq))
710ALLOCATE(tsoil_val(nsoil))
711#ifdef NEWPHYS
712ALLOCATE(isoil_val(nsoil))
713ALLOCATE(dsoil_val(nsoil))
714#endif
715ALLOCATE(wq2(ngrid,nlayer+1))    !!!!!
716ALLOCATE(wqsurf(ngrid,nq))       !!!!!
717ALLOCATE(wtsoil(ngrid,nsoil))    !!!!!
718#ifdef NEWPHYS
719ALLOCATE(wfluxrad(ngrid))
720ALLOCATE(wwstar(ngrid))
721ALLOCATE(wisoil(ngrid,nsoil))    !!!!!
722ALLOCATE(wdsoil(ngrid,nsoil))    !!!!!
723#endif
724ALLOCATE(pplev(ngrid,nlayer+1))  !!!!!
725ALLOCATE(pplay(ngrid,nlayer))    !!!!!
726ALLOCATE(pphi(ngrid,nlayer))     !!!!!
727ALLOCATE(pu(ngrid,nlayer))       !!!!!
728ALLOCATE(pv(ngrid,nlayer))       !!!!!
729ALLOCATE(pt(ngrid,nlayer))       !!!!!
730ALLOCATE(pw(ngrid,nlayer))       !!!!!
731ALLOCATE(pq(ngrid,nlayer,nq))    !!!!!
732! interm
733ALLOCATE(dz8w_prof(nlayer))
734ALLOCATE(p8w_prof(nlayer))
735ALLOCATE(p_prof(nlayer))
736ALLOCATE(t_prof(nlayer))
737ALLOCATE(t8w_prof(nlayer))
738ALLOCATE(u_prof(nlayer))
739ALLOCATE(v_prof(nlayer))
740ALLOCATE(z_prof(nlayer))
741!ALLOCATE(th_prof(nlayer))
742!ALLOCATE(rho_prof(nlayer))
743!ALLOCATE(pi_prof(nlayer))
744#ifdef NEWPHYS
745ALLOCATE(q_prof(nlayer,nq))
746ALLOCATE(wtnom(nq))              !!!!!
747#else
748ALLOCATE(water_vapor_prof(nlayer))
749ALLOCATE(water_ice_prof(nlayer))
750#endif
751
752
753!!!!!!!!!!!!!!!!!!!!!!!!!!!!
754!!!!!!!!!!!!!!!!!!!!!!!!!!!!
755!! PREPARE PHYSICS INPUTS !! 
756!!!!!!!!!!!!!!!!!!!!!!!!!!!!
757!!!!!!!!!!!!!!!!!!!!!!!!!!!!
758
759!SELECT CASE (MARS_MODE) !! ONLY ALLOW FOR MODES DEFINED IN Registry.EM
760!   CASE(4-10,13-19,22-41,43:)      !! -- CHANGE THIS if YOU ADDED CASES in REGISTRY.EM
761!   PRINT *, 'NOT SUPPORTED, to be done'
762!   STOP
763!END SELECT
764!!!!!!!!!!!!!!!!!!! FOR REFERENCE ; FROM REGISTRY.EM
765!package   default      mars==0                      -              -
766!package   water        mars==1                      -              scalar:qh2o,qh2o_ice
767!package   dust1        mars==2                      -              scalar:qdust
768!package   dust2eq      mars==3                      -              scalar:qdust,qdustn
769!package   newwater     mars==11                     -              scalar:qh2o,qh2o_ice,qdust,qdustn
770!package   radioac      mars==20                     -              scalar:qtrac1
771!package   radioac2     mars==21                     -              scalar:upward,downward
772!package   photochem    mars==42                     -              scalar:qco2,chem_co,chem_o,chem_o1d,chem_o2,chem_o3,chem_h,chem_h2,chem_oh,chem_ho2,chem_h2o2,chem_ch4,chem_n2,chem_ar,qh2o_ice,qh2o,qdust,qdustn
773!!!!!!!!!!!!!!!!!!! FOR REFERENCE
774
775#ifdef NEWPHYS
776!!! name of tracers to mimic entries in tracer.def
777!!! ----> IT IS IMPORTANT TO KEEP SAME ORDERING AS IN THE REGISTRY !!!!
778!!!                               SAME NAMING   AS IN THE LMD PHYSICS !!!!
779SELECT CASE (MARS_MODE)
780    CASE(0,10) 
781      wtnom(nq) = 'co2'
782    CASE(1)
783      wtnom(1)  = 'h2o_vap'
784      wtnom(2)  = 'h2o_ice'
785    CASE(2)
786      wtnom(1)  = 'dust01'     
787    CASE(3)
788      wtnom(1)  = 'dust_mass'
789      wtnom(2)  = 'dust_number'
790    CASE(11)
791      wtnom(1)  = 'h2o_vap'
792      wtnom(2)  = 'h2o_ice'
793      wtnom(3)  = 'dust_mass'
794      wtnom(4)  = 'dust_number'
795    CASE(12)
796      wtnom(1)  = 'h2o_vap'
797      wtnom(2)  = 'h2o_ice'
798      wtnom(3)  = 'dust_mass'
799      wtnom(4)  = 'dust_number'
800      wtnom(5)  = 'ccn_mass'
801      wtnom(6)  = 'ccn_number'
802    CASE(20)
803      wtnom(1) = 'qtrac1'
804    CASE(21)
805      wtnom(1) = 'upward'
806      wtnom(2) = 'downward'
807    CASE(42)
808      wtnom(1)  = 'co2'
809      wtnom(2)  = 'co'
810      wtnom(3)  = 'o'
811      wtnom(4)  = 'o1d'
812      wtnom(5)  = 'o2'
813      wtnom(6)  = 'o3'
814      wtnom(7)  = 'h'
815      wtnom(8)  = 'h2'
816      wtnom(9)  = 'oh'
817      wtnom(10)  = 'ho2'
818      wtnom(11)  = 'h2o2'
819      wtnom(12)  = 'ch4'
820      wtnom(13)  = 'n2'
821      wtnom(14)  = 'ar'
822      wtnom(15)  = 'h2o_ice'
823      wtnom(16)  = 'h2o_vap'
824      wtnom(17)  = 'dust_mass'
825      wtnom(18)  = 'dust_number'
826END SELECT
827#endif
828
829!!*******************************************!!
830!!*******************************************!!
831
832DO j = jps,jpe
833DO i = ips,ipe
834
835!!*******************************************!!
836!! FROM 3D WRF FIELDS TO 1D PHYSICS PROFILES !!
837!!*******************************************!!
838
839!-----------------------------------!
840! 1D subscript for physics "cursor" !
841!-----------------------------------!
842subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
843
844!--------------------------------------!
845! 1D columns : inputs for the physics  !
846! ... vert. coord., meteor. fields ... !
847!--------------------------------------!
848dz8w_prof(:) = dz8w(i,kps:kpe,j)   ! dz between full levels (m)   
849p8w_prof(:) = p8w(i,kps:kpe,j)     ! pressure full level (Pa) >> pplev
850p_prof(:) = p(i,kps:kpe,j)         ! pressure half level (Pa) >> pplay
851t_prof(:) = t(i,kps:kpe,j)         ! temperature half level (K) >> pt
852t8w_prof(:) = t8w(i,kps:kpe,j)     ! temperature full level (K)
853u_prof(:) = u(i,kps:kpe,j)         ! zonal wind (A-grid: unstaggered) half level (m/s) >> pu 
854v_prof(:) = v(i,kps:kpe,j)         ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv
855z_prof(:) = z(i,kps:kpe,j)         ! geopotential height half level (m) >> pphi/g
856!pi_prof(:) = exner(i,kps:kpe,j)    ! exner function (dimensionless) half level
857!rho_prof(:) = rho(i,kps:kpe,j)     ! density half level 
858!th_prof(:) = th(i,kps:kpe,j)       ! pot. temperature half level (K)
859
860!--------------------------------!
861! specific treatment for tracers !
862!--------------------------------!
863#ifndef NOPHYS
864#ifdef NEWPHYS
865IF (MARS_MODE .EQ. 0) THEN
866    q_prof(:,1)=0.95
867ELSE
868    q_prof(:,1:nq) = SCALAR(i,kps:kpe,j,2:nq+1)  !! the names were set above !! one dummy tracer in WRF
869ENDIF
870  !!!! CAS DU CO2
871  !DO iii=1,nq
872  ! IF ( wtnom(iii) .eq. 'co2' .and. (.not. restart)) q_prof(:,iii) = 0.95
873  !ENDDO
874  IF ((MARS_MODE .EQ. 20) .OR. (MARS_MODE .EQ. 21)) THEN
875   IF (firstcall .EQV. .true. .and. (.not. restart)) THEN
876      q_prof(:,:) = 0.95
877   ENDIF
878  ENDIF
879#else
880SELECT CASE (MARS_MODE)
881    CASE(0)  !! NO TRACERS (mars=0)
882    water_vapor_prof(:) = 0.
883    water_ice_prof(:) = 0.
884    CASE(1)  !! WATER CYCLE (mars=1)
885    water_vapor_prof(:) = scalar(i,kps:kpe,j,2)  !! H2O vapor is tracer 1 in the Registry for mars=1
886    water_ice_prof(:) = scalar(i,kps:kpe,j,3)    !! H2O ice is tracer 2 in the Registry for mars=1
887    CASE(2)  !! DUST CYCLE (mars=2)
888    water_vapor_prof(:) = 0.
889    water_ice_prof(:) = 0.
890    CASE(3:)
891    print *, 'OOOPS... not ready yet.'
892    STOP
893END SELECT
894#endif
895#endif
896
897!!**********************************************************!!
898!! STATIC FIELDS FOR THE PHYSICS - NEEDED ONLY AT FIRSTCALL !!
899!!**********************************************************!!
900needed_ini_phys : IF (firstcall .EQV. .true.) THEN
901
902!----------------------------------------!
903! Surface of each part of the grid (m^2) !
904!----------------------------------------!
905!aire_val = dx*dy                           !! 1. idealized cases - computational grid
906aire_val = (dx/msft(i,j))*(dy/msft(i,j))    !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
907!aire_val=dx*dy/msfu(i,j)                   !! 3. special for Mercator GCM-like simulations
908
909!---------------------------------------------!
910! Mass-point latitude and longitude (radians) !
911!---------------------------------------------!
912IF (JULYR .ne. 9999) THEN
913 lat_val = XLAT(i,j)*DEGRAD
914 lon_val = XLONG(i,j)*DEGRAD
915ELSE
916 !!! IDEALIZED CASE
917 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lat: ',lat_input
918 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lon: ',lon_input
919 lat_val = lat_input*DEGRAD
920 lon_val = lon_input*DEGRAD
921ENDIF
922
923!-----------------------------------------!
924! Gravity wave parametrization            !
925! NB: usually 0 in mesoscale applications !
926!-----------------------------------------!
927IF (JULYR .ne. 9999) THEN
928 zmea_val=M_GW(i,1,j)
929 zstd_val=M_GW(i,2,j)
930 zsig_val=M_GW(i,3,j)
931 zgam_val=M_GW(i,4,j)
932 zthe_val=M_GW(i,5,j)
933ELSE
934 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION GWdrag OFF'
935 zmea_val=0.
936 zstd_val=0.
937 zsig_val=0.
938 zgam_val=0.
939 zthe_val=0.
940ENDIF
941
942!---------------------------------!
943! Ground albedo & Thermal Inertia !
944!---------------------------------!
945IF (JULYR .ne. 9999) THEN
946 IF (CST_AL == 0) THEN
947 albedodat_val=M_ALBEDO(i,j)
948 ELSE
949 albedodat_val=CST_AL
950 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT ALBEDO ', albedodat_val
951 ENDIF
952 IF (CST_TI == 0) THEN
953 inertiedat_val=M_TI(i,j)
954 ELSE
955 inertiedat_val=CST_TI
956 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT THERMAL INERTIA ', inertiedat_val
957 ENDIF
958ELSE
959 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION albedo: ', CST_AL
960 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION inertia: ',CST_TI
961 albedodat_val=CST_AL
962 inertiedat_val=CST_TI
963ENDIF
964
965#ifdef NEWPHYS
966!----------------------------!
967! Variable surface roughness !
968!----------------------------!
969IF (JULYR .ne. 9999) THEN
970 IF (CST_Z0 == 0) THEN
971   z0_val = M_Z0(i,j)
972 ELSE
973   z0_val = CST_Z0
974   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT SURF ROUGHNESS (m) ',CST_Z0
975 ENDIF
976ELSE
977 IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION z0 (m) ', CST_Z0
978 z0_val=CST_Z0
979ENDIF
980!!!! ADDITIONAL SECURITY. THIS MIGHT HAPPEN WITH OLD INIT FILES.
981IF (z0_val == 0.) THEN
982   IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
983   z0_val = 0.01
984ENDIF
985!!!! ADDITIONAL SECURITY. INTERP+SMOOTH IN GEOGRID MIGHT YIELD NEGATIVE Z0 !!!
986IF (z0_val < 0.) THEN
987   PRINT *, 'WELL, z0 is NEGATIVE, this is impossible. better stop here.'
988   PRINT *, 'advice --> correct interpolation / smoothing of z0 in WPS'
989   PRINT *, '           -- or check the constant value set in namelist.input'
990   STOP
991ENDIF
992#endif
993
994!---------------------------------------------------------!
995! Ground geopotential                                     !
996! (=g*HT(i,j), only used in the microphysics: newcondens) !
997!---------------------------------------------------------!
998phisfi_val=g*(z_prof(1)-dz8w_prof(1)/2.)   !! NB: z_prof are half levels, so z_prof(1) is not the surface
999
1000!-----------------------------------------------!
1001! Ground temperature, emissivity, CO2 ice cover !
1002!-----------------------------------------------!
1003IF (M_TSURF(i,j) .gt. 0.) THEN
1004  tsurf_val=M_TSURF(i,j)
1005ELSE
1006  tsurf_val=TSK(i,j) ! retro-compatibility
1007ENDIF
1008emis_val=M_EMISS(i,j)
1009co2ice_val=M_CO2ICE(i,j)
1010
1011!------------------------!
1012! Deep soil temperatures !
1013!------------------------!
1014IF (JULYR .ne. 9999) THEN
1015  IF (M_TSOIL(i,1,j) .gt. 0.) THEN
1016   tsoil_val(:)=M_TSOIL(i,:,j)
1017  ELSE
1018   tsoil_val = tsoil_val*0. + tsurf_val
1019  ENDIF
1020#ifdef NEWPHYS
1021  isoil_val(:)=M_ISOIL(i,:,j)
1022  dsoil_val(:)=M_DSOIL(i,:,j)
1023#endif
1024ELSE
1025  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION tsoil is set to tsurf'
1026  do k=1,nsoil
1027   IF (.not.restart) THEN
1028     tsoil_val(k) = tsurf_val
1029   ELSE
1030     !this is a restart run. We must not set tsoil to tsurf in the init.
1031     !tsoil was saved in physiq.F under the name M_TSOIL in the restart file
1032     !(see Registry)
1033     tsoil_val(k)=M_TSOIL(i,k,j)
1034   ENDIF
1035
1036#ifdef NEWPHYS
1037   IF ( nsoil .lt. 18 ) THEN
1038       PRINT *,'** Mars ** WRONG NUMBER OF SOIL LAYERS. SHOULD BE 18 AND IT IS ',nsoil
1039       STOP
1040   ENDIF
1041   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION isoil and dsoil standard'
1042   isoil_val(k) = inertiedat_val
1043!  dsoil_val(k) = sqrt(887.75/3.14)*((2.**(k-0.5))-1.) * inertiedat_val / wvolcapa    ! old setting
1044   dsoil_val(k) = 2.E-4 * (2.**(k-0.5-1.))                                            ! new gcm settings
1045   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** ISOIL DSOIL are ',isoil_val(k), dsoil_val(k)
1046#endif
1047  enddo
1048ENDIF
1049
1050!-------------------!
1051! Tracer at surface !
1052!-------------------!
1053#ifndef NOPHYS
1054SELECT CASE (MARS_MODE)
1055    CASE(0) 
1056    qsurf_val(:)=0.
1057    CASE(1)
1058    qsurf_val(1)=0.
1059    qsurf_val(2)=M_H2OICE(i,j)    !! logique avec wtnom(2) = 'h2o_ice' defini ci-dessus
1060                                   !! ----- retrocompatible ancienne physique
1061                                   !! ----- [H2O ice is last tracer in qsurf in LMD physics]
1062    CASE(2) 
1063    qsurf_val(1)=0.                !! not coupled with lifting for the moment [non remobilise]
1064#ifdef NEWPHYS
1065    CASE(3)
1066    qsurf_val(1)=q_prof(1,1)                !!! temporaire, a definir       
1067    qsurf_val(2)=q_prof(1,2)     
1068    CASE(11)
1069    qsurf_val(1)=0.
1070    qsurf_val(2)=M_H2OICE(i,j)    !! logique avec wtnom(2) = 'h2o_ice' defini ci-dessus
1071    qsurf_val(3)=0.                !! not coupled with lifting for the moment [non remobilise]
1072    qsurf_val(4)=0.
1073    CASE(12)   
1074    qsurf_val(1)=0.
1075    qsurf_val(2)=M_H2OICE(i,j)    !! logique avec wtnom(2) = 'h2o_ice' defini ci-dessus
1076    qsurf_val(3)=0.                !! not coupled with lifting for the moment [non remobilise]
1077    qsurf_val(4)=0.
1078    qsurf_val(5)=0.
1079    qsurf_val(6)=0.
1080    CASE(20)
1081    qsurf_val(:)=0.
1082    CASE(21)
1083    qsurf_val(:)=0.
1084#else
1085    CASE(3:)
1086    print *, 'OOOPS... not ready yet.'
1087    STOP
1088#endif
1089
1090END SELECT
1091#endif
1092
1093!-------------------!
1094! Slope inclination !
1095!-------------------!
1096IF (JULYR .ne. 9999) THEN
1097  theta_val=atan(sqrt( (1000.*SLPX(i,j))**2 + (1000.*SLPY(i,j))**2 ))
1098  theta_val=theta_val/DEGRAD
1099ELSE
1100  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION slope inclination is 0'
1101  theta_val=0.
1102ENDIF
1103
1104!-------------------------------------------!
1105! Slope orientation; 0 is north, 90 is east !
1106!-------------------------------------------!
1107IF (JULYR .ne. 9999) THEN
1108  psi_val=-90.*DEGRAD-atan(SLPY(i,j)/SLPX(i,j))
1109  if (SLPX(i,j) .ge. 0.) then
1110    psi_val=psi_val-180.*DEGRAD
1111  endif
1112  psi_val=360.*DEGRAD+psi_val
1113  psi_val=psi_val/DEGRAD
1114  psi_val = MODULO(psi_val+180.,360.)
1115ELSE
1116  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION slope orientation is 0 (well, whatever)'
1117  psi_val=0.
1118ENDIF
1119
1120!-----------------!
1121! Optional output !
1122!-----------------!
1123IF ( (i == ips) .AND. (j == jps) ) THEN
1124PRINT *,'lat/lon ', lat_val/DEGRAD, lon_val/DEGRAD
1125PRINT *,'emiss ', emis_val
1126PRINT *,'albedo ', albedodat_val
1127PRINT *,'inertie ', inertiedat_val
1128PRINT *,'phi ',phisfi_val
1129PRINT *,'tsurf ',tsurf_val
1130PRINT *,'aire ',aire_val
1131PRINT *,'z_prof ',z_prof
1132PRINT *,'dz8w_prof',dz8w_prof
1133PRINT *,'p8w_prof ',p8w_prof
1134PRINT *,'p_prof ',p_prof
1135PRINT *,'t_prof ',t_prof
1136PRINT *,'t8w_prof ',t8w_prof
1137PRINT *,'u_prof ',u_prof
1138PRINT *,'v_prof ',v_prof
1139PRINT *,'tsoil ',tsoil_val
1140PRINT *,'qsurf ',qsurf_val
1141#ifdef NEWPHYS
1142PRINT *,'isoil ',isoil_val
1143PRINT *,'dsoil ',dsoil_val
1144PRINT *,'q_prof ',q_prof
1145PRINT *,'z0 ',z0_val
1146#endif
1147ENDIF
1148
1149!-------------------------!
1150!-------------------------!
1151!      PROVISOIRE         !
1152!-------------------------!
1153!-------------------------!
1154IF (.not. restart) THEN
1155   q2_val(:)=1.E-6      !PBL wind variance
1156#ifdef NEWPHYS
1157   fluxrad_val=0.
1158   wstar_val=0.
1159#endif
1160ELSE
1161   q2_val(:)=M_Q2(i,:,j)
1162#ifdef NEWPHYS
1163   fluxrad_val=M_FLUXRAD(i,j)
1164   wstar_val=M_WSTAR(i,j)
1165#endif
1166ENDIF
1167
1168!-----------------!
1169! Fill the arrays !
1170!-----------------!
1171aire_vec(subs) = aire_val   !! NB: total area in square km is SUM(aire_vec)/1.0E6
1172lat_vec(subs) = lat_val
1173lon_vec(subs) = lon_val
1174wphisfi(subs) = phisfi_val
1175walbedodat(subs) = albedodat_val
1176winertiedat(subs) = inertiedat_val
1177wzmea(subs) = zmea_val
1178wzstd(subs) = zstd_val
1179wzsig(subs) = zsig_val
1180wzgam(subs) = zgam_val
1181wzthe(subs) = zthe_val
1182wtsurf(subs) = tsurf_val
1183wco2ice(subs) = co2ice_val
1184wemis(subs) = emis_val
1185wq2(subs,:) = q2_val(:)
1186wqsurf(subs,:) = qsurf_val(:)
1187wtsoil(subs,:) = tsoil_val(:)
1188#ifdef NEWPHYS
1189wfluxrad(subs) = fluxrad_val
1190wwstar(subs) = wstar_val
1191wisoil(subs,:) = isoil_val(:)
1192wdsoil(subs,:) = dsoil_val(:)
1193wz0tab(subs) = z0_val
1194#endif
1195wtheta(subs) = theta_val
1196wpsi(subs) = psi_val
1197
1198ENDIF needed_ini_phys
1199
1200!!*****************************!!
1201!! PREPARE "FOOD" FOR PHYSIQ.F !!
1202!!*****************************!!
1203
1204!---------------------------------------------!
1205! in LMD physics, geopotential must be        !
1206! expressed with respect to the local surface !
1207!---------------------------------------------!
1208pphi(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) )
1209
1210!--------------------------------!
1211! Dynamic fields for LMD physics !
1212!--------------------------------!
1213pplev(subs,1:nlayer) = p8w_prof(1:nlayer)  !! NB: last level: no data
1214pplay(subs,:) = p_prof(:)
1215pt(subs,:) = t_prof(:)
1216pu(subs,:) = u_prof(:)
1217pv(subs,:) = v_prof(:)
1218pw(subs,:) = 0   !! NB: not used in the physics, only diagnostic...
1219!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1220!! for IDEALIZED CASES ONLY
1221IF (JULYR .eq. 9999) pplev(subs,nlayer+1)=0.  !! pplev(subs,nlayer+1)=ptop >> NO !
1222!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1223
1224! NOTE:
1225! IF ( pplev(subs,nlayer+1) .le. 0 ) pplev(subs,nlayer+1)=ptop
1226! cree des diagnostics delirants et aleatoires dans le transfert radiatif
1227
1228!---------!
1229! Tracers ! 
1230!---------!
1231#ifndef NOPHYS
1232#ifdef NEWPHYS
1233pq(subs,:,:) = q_prof(:,:)  !! traceurs generiques, seuls noms sont specifiques
1234#else
1235SELECT CASE (MARS_MODE)
1236    CASE(0)  !! NO TRACERS (mars=0)
1237    pq(subs,:,nq) = water_vapor_prof(:) !! NB: which is 0, actually (see above)
1238    CASE(1)  !! WATER CYCLE (mars=1)
1239    pq(subs,:,nq) = water_vapor_prof(:)
1240    pq(subs,:,nq-1) = water_ice_prof(:)
1241    CASE(2)  !! DUST CYCLE (mars=2)
1242    pq(subs,:,nq) = water_vapor_prof(:) !! NB: which is 0, actually (see above)
1243    CASE(3:)
1244    print *, 'OOOPS... not ready yet.'
1245    STOP   
1246END SELECT
1247#endif
1248#endif
1249ENDDO
1250ENDDO
1251
1252!!*****************!!
1253!! CLEAN THE PLACE !!
1254!!*****************!!
1255DEALLOCATE(q2_val)
1256DEALLOCATE(qsurf_val)
1257DEALLOCATE(tsoil_val)
1258#ifdef NEWPHYS
1259DEALLOCATE(isoil_val)
1260DEALLOCATE(dsoil_val)
1261#endif
1262DEALLOCATE(dz8w_prof)
1263DEALLOCATE(z_prof)
1264DEALLOCATE(p8w_prof)
1265DEALLOCATE(p_prof)
1266DEALLOCATE(t_prof)
1267DEALLOCATE(t8w_prof)
1268DEALLOCATE(u_prof)
1269DEALLOCATE(v_prof)
1270#ifdef NEWPHYS
1271DEALLOCATE(q_prof)
1272#else
1273DEALLOCATE(water_vapor_prof)
1274DEALLOCATE(water_ice_prof)
1275#endif
1276    !!! no use
1277    !DEALLOCATE(pi_prof)
1278    !DEALLOCATE(rho_prof)
1279    !DEALLOCATE(th_prof)
1280
1281
1282!!!!!!!!!!!!!!!!!!!!!!
1283!!!!!!!!!!!!!!!!!!!!!!
1284!! CALL LMD PHYSICS !! 
1285!!!!!!!!!!!!!!!!!!!!!!
1286!!!!!!!!!!!!!!!!!!!!!!
1287
1288
1289!!***********************!!
1290!! INIFIS, AT FIRST CALL !!
1291!!***********************!!
1292IF (firstcall .EQV. .true.) THEN
1293#ifndef NOPHYS
1294print *, '** Mars ** LMD INITIALIZATION'
1295#include "../call_meso_inifis.inc"
1296!!! le # est important pour newphys
1297#endif
1298DEALLOCATE(aire_vec)
1299DEALLOCATE(lat_vec)
1300DEALLOCATE(lon_vec)
1301DEALLOCATE(walbedodat)
1302DEALLOCATE(winertiedat)
1303DEALLOCATE(wphisfi)
1304DEALLOCATE(wzmea)
1305DEALLOCATE(wzstd)
1306DEALLOCATE(wzsig)
1307DEALLOCATE(wzgam)
1308DEALLOCATE(wzthe)
1309DEALLOCATE(wtheta)
1310DEALLOCATE(wpsi)
1311#ifdef NEWPHYS
1312DEALLOCATE(wz0tab)
1313#endif
1314ENDIF
1315
1316!!********!!
1317!! PHYSIQ !!
1318!!********!!
1319call_physics : IF (wappel_phys .ne. 0.) THEN
1320
1321!-------------------------------------------------------------------------------!
1322! outputs:                                                                      !       
1323!    pdu(ngrid,nlayermx)        \                                               !
1324!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding    !
1325!    pdt(ngrid,nlayermx)         /  variables due to physical processes.        !
1326!    pdq(ngrid,nlayermx)        /                                               !
1327!    pdpsrf(ngrid)             /                                                !
1328!    tracerdyn                 call tracer in dynamical part of GCM ?           !
1329!-------------------------------------------------------------------------------!
1330pdpsrf(:)=0.
1331pdu(:,:)=0.
1332pdv(:,:)=0.
1333pdt(:,:)=0.
1334pdq(:,:,:)=0.
1335#ifndef NOPHYS
1336print *, '** Mars ** CALL TO LMD PHYSICS'
1337#include "../call_meso_physiq.inc"
1338!!! le # est important pour newphys
1339#endif
1340DEALLOCATE(pplev)
1341DEALLOCATE(pplay)
1342DEALLOCATE(pphi)
1343DEALLOCATE(pu)
1344DEALLOCATE(pv)
1345DEALLOCATE(pt)
1346DEALLOCATE(pw)
1347DEALLOCATE(pq)
1348!DEALLOCATE(wtsurf)
1349!DEALLOCATE(wco2ice)
1350DEALLOCATE(wemis)
1351!DEALLOCATE(wq2)
1352!DEALLOCATE(wqsurf)
1353!DEALLOCATE(wtsoil)
1354!DEALLOCATE(wwstar)
1355!DEALLOCATE(wfluxrad)
1356#ifdef NEWPHYS
1357DEALLOCATE(wisoil)
1358DEALLOCATE(wdsoil)
1359DEALLOCATE(wtnom)
1360#endif
1361
1362
1363!-------------------------------!
1364! PHYSIQ OUTPUT IN THE WRF FILE !
1365!-------------------------------!
1366#ifndef NOPHYS
1367DO j = jps,jpe
1368DO i = ips,ipe
1369subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
1370
1371!!!! SEMBLE T IL PROBLEME AVEC NEWPHYS .... qui marche tres bien sans !
1372!#ifndef NEWPHYS           
1373!            !!! sparadrap pour regler un probleme avec mpiifort en LES
1374!            !!! -- HFX apparaissait nul aux interfaces des tuiles
1375!            if ( output_tab2d(subs,ind_HFX) .le. 1.e-8 ) then
1376!               !print *, 'HFX is zero !!! ', i, j
1377!               !print *, 'I substituted the value right next to it ', output_tab2d(subs+1,ind_HFX)
1378!               output_tab2d(subs,ind_HFX) = output_tab2d(subs+1,ind_HFX)
1379!            endif
1380!#endif
1381
1382#include "module_lmd_driver_output3.inc" 
1383       !  ^-- generated from Registry
1384ENDDO
1385ENDDO
1386DEALLOCATE(output_tab2d)
1387DEALLOCATE(output_tab3d)
1388#endif
1389
1390!!!!!! PATCH SPECIAL STORM
1391#ifdef NEWPHYS
1392#ifdef storm
1393#include "../mars_lmd/storm.inc"
1394#endif
1395#endif
1396!!!!!! PATCH SPECIAL STORM
1397
1398
1399!---------------------------------------------------------------------------------!
1400! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
1401!---------------------------------------------------------------------------------!
1402#ifdef SPECIAL_NEST_SAVE
1403dp_save(:,id)=pdpsrf(:)
1404du_save(:,:,id)=pdu(:,:)
1405dv_save(:,:,id)=pdv(:,:)
1406dt_save(:,:,id)=pdt(:,:)
1407dq_save(:,:,:,id)=pdq(:,:,:)
1408#ifndef NORESTART
1409save_tsoil_restart(:,:,id)=wtsoil(:,:)
1410save_co2ice_restart(:,id)=wco2ice(:)
1411save_q2_restart(:,:,id)=wq2(:,:)
1412save_qsurf_restart(:,:,id)=wqsurf(:,:)
1413save_tsurf_restart(:,id)=wtsurf(:)
1414#ifdef NEWPHYS
1415save_wstar_restart(:,id)=wwstar(:)
1416save_fluxrad_restart(:,id)=wfluxrad(:)
1417#endif
1418#endif
1419#else
1420dp_save(:)=pdpsrf(:)
1421du_save(:,:)=pdu(:,:)
1422dv_save(:,:)=pdv(:,:)
1423dt_save(:,:)=pdt(:,:)
1424dq_save(:,:,:)=pdq(:,:,:)
1425#ifndef NORESTART
1426save_tsoil_restart(:,:)=wtsoil(:,:)
1427save_co2ice_restart(:)=wco2ice(:)
1428save_q2_restart(:,:)=wq2(:,:)
1429save_qsurf_restart(:,:)=wqsurf(:,:)
1430save_tsurf_restart(:)=wtsurf(:)
1431#ifdef NEWPHYS
1432save_wstar_restart(:)=wwstar(:)
1433save_fluxrad_restart(:)=wfluxrad(:)
1434#endif
1435#endif
1436#endif
1437DEALLOCATE(wtsoil)
1438DEALLOCATE(wco2ice)
1439DEALLOCATE(wq2)
1440DEALLOCATE(wqsurf)
1441DEALLOCATE(wtsurf)
1442#ifdef NEWPHYS
1443DEALLOCATE(wfluxrad)
1444DEALLOCATE(wwstar)
1445#endif
1446ENDIF call_physics
1447
1448ENDIF    ! end of BIG LOOP 2.
1449!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1451!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1452!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1453
1454
1455!!***************************!!
1456!! DEDUCE TENDENCIES FOR WRF !!
1457!!***************************!!
1458RTHBLTEN(ims:ime,kms:kme,jms:jme)=0.
1459RUBLTEN(ims:ime,kms:kme,jms:jme)=0.
1460RVBLTEN(ims:ime,kms:kme,jms:jme)=0.
1461PSFC(ims:ime,jms:jme)=p8w(ims:ime,kms,jms:jme) ! was done in surface driver in regular WRF
1462!------------------------------------------------------------------!
1463! WRF adds the physical and dynamical tendencies                   !
1464! thus the tendencies must be passed as 'per second' tendencies    !
1465! --when physics is not called, the applied physical tendency      !
1466! --is the one calculated during the last call to physics          !
1467!------------------------------------------------------------------!
1468
1469! PBL top index reading for MARS_MODE 21 :
1470IF (MARS_MODE .EQ. 21) THEN
1471      IF (firstcall .EQV. .true.) THEN
1472           open(unit=15,file='input_zipbl',form='formatted',status='old')
1473           rewind(15)
1474           end_of_file = .false.
1475           kk = 0
1476           do while (.not. end_of_file)
1477             read(15,*,end=101) zipbl(kk+1),zilct(kk+1)
1478             write(*,*) kk, zipbl(kk+1),zilct(kk+1)
1479             kk = kk+1
1480             go to 112
1481101          end_of_file = .true.
1482112          continue
1483           enddo
1484           close(unit=15,status = 'keep')
1485      ENDIF
1486      lct=lct_input + elaps/3700.
1487      zilctmin=MINVAL(ABS(zilct(:)-lct))
1488      lctindex=0
1489      DO kk=1,500
1490        IF(ABS(zilct(kk)-lct) .eq. zilctmin) lctindex=kk
1491      ENDDO
1492      IF(lctindex .eq. 0) print*, 'probleme index'
1493        IF ((zilct(lctindex) .gt. 12.) .and.( zilct(lctindex) .lt. 18)) THEN
1494            ziindex=MAX(0.,zipbl(lctindex)-7.)
1495        ELSE
1496            ziindex=0.
1497        ENDIF
1498ENDIF
1499
1500#ifdef NOPHYS
1501!!!!!!!!!!!!!!!!!!!
1502!!!!!!!!!!!!!!!!!!!
1503IF (firstcall .EQV. .true.) THEN
1504  ALLOCATE(dtrad(nlayer+1))
1505  DO k=kps,kpe+1
1506    !! the only solution to avoid 0 points in M_Q2
1507    !! -- in each domain decomposition case
1508    dtrad(k) = MAXVAL(M_Q2(:,k,:)) + MINVAL(M_Q2(:,k,:))
1509  ENDDO
1510  print*,"HEATING RATE",dtrad(kps:kpe+1)
1511ENDIF
1512DO i= 1,ngrid
1513  pdt(i,kps:kpe) = dtrad(kps:kpe)
1514ENDDO
1515!!!!!!!!!!!!!!!!!!!
1516!!!!!!!!!!!!!!!!!!!
1517#endif
1518
1519DO j = jps,jpe
1520DO i = ips,ipe
1521subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
1522
1523!------------!
1524! zonal wind !
1525!------------!
1526RUBLTEN(i,kps:kpe,j) = pdu(subs,kps:kpe)                       
1527
1528!-----------------!
1529! meridional wind !
1530!-----------------!
1531RVBLTEN(i,kps:kpe,j) = pdv(subs,kps:kpe)                       
1532
1533!-----------------------!
1534! potential temperature !
1535!-----------------------!
1536! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
1537RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe) / exner(i,kps:kpe,j)   
1538
1539!---------------------------!
1540! update surface pressure   !
1541! (cf CO2 cycle in physics) !
1542!---------------------------!
1543PSFC(i,j)=PSFC(i,j)+pdpsrf(subs)*dt    !!! here dt is needed
1544
1545!------------------------------------!
1546! Save key variables for restart ! 
1547!------------------------------------!
1548#ifndef NORESTART
1549#ifdef SPECIAL_NEST_SAVE
1550M_TSOIL(i,:,j)=save_tsoil_restart(subs,:,id)
1551M_CO2ICE(i,j)=save_co2ice_restart(subs,id)
1552M_Q2(i,kps:kpe+1,j)=save_q2_restart(subs,:,id)
1553SELECT CASE (MARS_MODE)
1554   CASE (1,11,12)
1555     M_H2OICE(i,j)=save_qsurf_restart(subs,2,id)  !! see above Tracer at surface
1556END SELECT
1557M_TSURF(i,j)=save_tsurf_restart(subs,id)
1558#ifdef NEWPHYS
1559M_WSTAR(i,j)=save_wstar_restart(subs,id)
1560M_FLUXRAD(i,j)=save_fluxrad_restart(subs,id)
1561#endif
1562#else
1563M_TSOIL(i,:,j)=save_tsoil_restart(subs,:)
1564M_CO2ICE(i,j)=save_co2ice_restart(subs)
1565M_Q2(i,kps:kpe+1,j)=save_q2_restart(subs,:)
1566SELECT CASE (MARS_MODE)
1567   CASE (1,11,12)
1568     M_H2OICE(i,j)=save_qsurf_restart(subs,2)  !! see above Tracer at surface
1569END SELECT
1570M_TSURF(i,j)=save_tsurf_restart(subs)
1571#ifdef NEWPHYS
1572M_WSTAR(i,j)=save_wstar_restart(subs)
1573M_FLUXRAD(i,j)=save_fluxrad_restart(subs)
1574#endif
1575#endif
1576#endif
1577
1578!---------!
1579! Tracers ! 
1580!---------!
1581#ifndef NOPHYS
1582#ifdef NEWPHYS
1583SCALAR(i,kps:kpe,j,1)=0.
1584
1585SELECT CASE (MARS_MODE)
1586CASE(0)
1587        SCALAR(i,kps:kpe,j,:)=0.
1588CASE(20)
1589   !! Mars mode 20 : add a passive tracer with radioactive-like decay
1590   IF ( (i == ips) .AND. (j == jps) )   print *, 'RADIOACTIVE-LIKE TRACER WITH SOURCE AT SURFACE LAYER.'
1591   tau_decay=60.*10. !! why not make it a namelist argument?
1592   SCALAR(i,kps:kpe,j,2) = SCALAR(i,kps:kpe,j,2)*exp(-dt/tau_decay)
1593   SCALAR(i,1,j,2) = SCALAR(i,1,j,2) + 1. !! this tracer is emitted in the surface layer
1594CASE(21)
1595   !! Mars mode 21 : add a passive tracer with radioactive-like decay at surface and pbl top
1596   IF ( (i == ips) .AND. (j == jps) )   print *, 'RADIOACTIVE-LIKE TRACER WITH SOURCE AT SURFACE LAYER AND PBL TOP.'
1597   tau_decay=60.*10. !! why not make it a namelist argument?
1598   SCALAR(i,kps:kpe,j,2) = SCALAR(i,kps:kpe,j,2)*exp(-dt/tau_decay)
1599   SCALAR(i,kps:kpe,j,3) = SCALAR(i,kps:kpe,j,3)*exp(-dt/tau_decay)
1600   SCALAR(i,1,j,2) = SCALAR(i,1,j,2) + 1. !! this tracer is emitted in the surface layer
1601   IF ( (i == ips) .AND. (j == jps) )   print *, 'index of pbl emission: ',ziindex
1602   IF (ziindex .NE. 0) THEN
1603   SCALAR(i,ziindex,j,3) = SCALAR(i,ziindex,j,3) + 1. !! this tracer is emitted at pbl top
1604   ENDIF
1605CASE DEFAULT
1606   SCALAR(i,kps:kpe,j,2:nq+1)=SCALAR(i,kps:kpe,j,2:nq+1)+pdq(subs,kps:kpe,1:nq)*dt  !!! here dt is needed
1607END SELECT
1608#else
1609SELECT CASE (MARS_MODE)
1610CASE(0)
1611        SCALAR(i,kps:kpe,j,:)=0.
1612CASE(1)
1613        !!! Water vapor
1614        SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*dt  !!! here dt is needed
1615        !!! Water ice
1616        SCALAR(i,kps:kpe,j,3)=SCALAR(i,kps:kpe,j,3)+pdq(subs,kps:kpe,nq-1)*dt  !!! here dt is needed
1617CASE(2)
1618        !!! Dust
1619        SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*dt  !!! here dt is needed
1620CASE(3:)
1621        print *, 'OOOPS... not ready yet.'
1622        STOP   
1623END SELECT
1624#endif
1625#endif
1626        !!TODO: check if adding the whole tendency once, and set the
1627        !!TODO: following tendencies to 0 until physics is called again
1628        !!TODO: is a good strategy ?
1629        !RUBLTEN(i,kps:kpe,j) = pdu(subs,kps:kpe)*ptimestep/dt
1630        !RVBLTEN(i,kps:kpe,j) = pdv(subs,kps:kpe)*ptimestep/dt
1631        !RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe)*ptimestep/dt
1632        !RTHBLTEN(i,kps:kpe,j) = RTHBLTEN(i,kps:kpe,j)/exner(i,kps:kpe,j)
1633        !PSFC(i,j)=PSFC(i,j)+pdpsrf(subs)*ptimestep/dt     
1634        !SELECT CASE (MARS_MODE)
1635        !CASE(0)
1636        !SCALAR(i,kps:kpe,j,:)=0.
1637        !CASE(1)
1638        !!!! Water vapor
1639        !SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*ptimestep/dt
1640        !!!! Water ice
1641        !SCALAR(i,kps:kpe,j,3)=SCALAR(i,kps:kpe,j,3)+pdq(subs,kps:kpe,nq-1)*ptimestep/dt
1642        !CASE(2)
1643        !!!! Dust
1644        !SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*ptimestep/dt
1645        !END SELECT
1646
1647ENDDO
1648ENDDO
1649DEALLOCATE(pdpsrf)
1650DEALLOCATE(pdu)
1651DEALLOCATE(pdv)
1652DEALLOCATE(pdt)
1653DEALLOCATE(pdq)
1654
1655!!---------!
1656!! display !
1657!!---------!
1658!PRINT *, '** Mars ** Results from LMD physics'
1659!PRINT *, 'u non-zero tendencies'
1660!!PRINT *, 'max',MAXVAL(RUBLTEN, MASK=RUBLTEN/=0.),&
1661!!        ' at',MAXLOC(RUBLTEN, MASK=RUBLTEN/=0.)
1662!!PRINT *, 'min',MINVAL(RUBLTEN, MASK=RUBLTEN/=0.),&
1663!!        ' at',MINLOC(RUBLTEN, MASK=RUBLTEN/=0.)
1664!PRINT *, RUBLTEN(10,1,10), RUBLTEN(10,15,10)
1665!PRINT *, 'v non-zero tendencies'
1666!!PRINT *, 'max',MAXVAL(RVBLTEN, MASK=RVBLTEN/=0.),&
1667!!        ' at',MAXLOC(RVBLTEN, MASK=RVBLTEN/=0.)
1668!!PRINT *, 'min',MINVAL(RVBLTEN, MASK=RVBLTEN/=0.),&
1669!!        ' at',MINLOC(RVBLTEN, MASK=RVBLTEN/=0.)
1670!PRINT *, RVBLTEN(10,1,10), RVBLTEN(10,15,10)
1671!!! STOP IF CRASH
1672!IF (MAXVAL(RUBLTEN, MASK=RUBLTEN/=0.) == 0.) STOP
1673!IF (MAXVAL(RVBLTEN, MASK=RVBLTEN/=0.) == 0.) STOP
1674
1675!!*****!!
1676!! END !!
1677!!*****!!
1678print *,'** Mars ** END LMD PHYSICS'
1679!----------------------------------------------------------------!
1680! use debug (see solve_em) if you wanna display some message ... !
1681!----------------------------------------------------------------!
1682END SUBROUTINE lmd_driver
1683
1684!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1685      real function ls2sol(ls)
1686
1687!c  Returns solar longitude, Ls (in deg.), from day number (in sol),
1688!c  where sol=0=Ls=0 at the northern hemisphere spring equinox
1689
1690      implicit none
1691
1692!c  Arguments:
1693      real ls
1694
1695!c  Local:
1696      double precision xref,zx0,zteta,zz
1697!c      xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
1698      double precision year_day
1699      double precision peri_day,timeperi,e_elips
1700      double precision pi,degrad
1701      parameter (year_day=668.6d0) ! number of sols in a amartian year
1702!c      data peri_day /485.0/
1703      parameter (peri_day=485.35d0) ! date (in sols) of perihelion
1704!c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
1705      parameter (timeperi=1.90258341759902d0)
1706      parameter (e_elips=0.0934d0)  ! eccentricity of orbit
1707      parameter (pi=3.14159265358979d0)
1708      parameter (degrad=57.2957795130823d0)
1709
1710      if (abs(ls).lt.1.0e-5) then
1711         if (ls.ge.0.0) then
1712            ls2sol = 0.0
1713         else
1714            ls2sol = year_day
1715         end if
1716         return
1717      end if
1718
1719      zteta = ls/degrad + timeperi
1720      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
1721      xref = zx0-e_elips*dsin(zx0)
1722      zz = xref/(2.*pi)
1723      ls2sol = zz*year_day + peri_day
1724      if (ls2sol.lt.0.0) ls2sol = ls2sol + year_day
1725      if (ls2sol.ge.year_day) ls2sol = ls2sol - year_day
1726
1727      return
1728      end function ls2sol
1729!!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1730
1731END MODULE module_lmd_driver
Note: See TracBrowser for help on using the repository browser.