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

Last change on this file since 3094 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
RevLine 
[11]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
[28]12!           - additions for new soil model in physics - January 2010
[65]13!           - unified module_lmd_driver: old, new phys and LES - February 2011
[11]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, &
[2335]30        num_3d_m,MOIST,&
[11]31        MARS_MODE, &
[1236]32        M_ALBEDO,M_TI,M_CO2ICE,M_EMISS, &
33        M_H2OICE, &
34        M_TSOIL, &
35        M_Q2, &
36        M_TSURF, &
[678]37#ifdef NEWPHYS
[1236]38        M_FLUXRAD, &
39        M_WSTAR, &
40        M_ISOIL, &
41        M_DSOIL, &
42        M_Z0, &
[315]43        CST_Z0, &
[28]44#endif
[1236]45        M_GW, &
[11]46        NUM_SOIL_LAYERS, &
47        CST_AL, &
48        CST_TI, &
[34]49        isfflx, &
50        diff_opt, &
51        km_opt, &
[11]52        HISTORY_INTERVAL, &
[156]53#ifndef NOPHYS
[11]54#include "module_lmd_driver_output1.inc"
[156]55#endif
[674]56        SLPX,SLPY,RESTART)
[11]57! NB: module_lmd_driver_output1.inc : output arguments generated from Registry
58
[34]59
60
61
62
[11]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
[156]78#ifndef NOPHYS
[11]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"
[156]94#endif
[11]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
[34]110INTEGER, INTENT(IN  ) :: isfflx,diff_opt,km_opt
[11]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,  &
[1236]118     M_ALBEDO,M_TI,M_EMISS, &
[11]119     SLPX,SLPY
[674]120REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT   )  :: &
[1236]121     M_CO2ICE,M_H2OICE, &
122     M_TSURF
[11]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
[1388]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 ) :: &
[1236]129     M_Q2
[11]130!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
131INTEGER, INTENT(IN   ) :: HISTORY_INTERVAL
132!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[674]133REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(INOUT )  :: &
[1236]134     M_TSOIL
[28]135#ifdef NEWPHYS
[315]136REAL, INTENT(IN  ) :: CST_Z0
[28]137REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(IN   )  :: &
[1236]138     M_ISOIL, M_DSOIL         
[315]139REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  :: &
[1236]140     M_Z0
[678]141REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT   )  :: &
[1236]142     M_FLUXRAD,M_WSTAR
[28]143#endif
[11]144REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: &
[1236]145     M_GW
[11]146! 4D arrays
[2335]147INTEGER, INTENT(IN ) :: num_3d_s,num_3d_m
[11]148REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_s ), INTENT(INOUT ) :: &
149     scalar
[2335]150REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_m ), INTENT(INOUT ) :: &
151     moist
[674]152! Logical
153LOGICAL, INTENT(IN ) :: restart
[11]154
155!-------------------------------------------
156! OUTPUT VARIABLES
157!-------------------------------------------
158!
159! Generated from Registry
160!
161! default definitions :
162! 2D : TSK, PSFC
163! 3D : RTHBLTEN,RUBLTEN,RVBLTEN
[156]164#ifndef NOPHYS
[11]165#include "module_lmd_driver_output2.inc"
166   REAL, DIMENSION(:,:), ALLOCATABLE :: output_tab2d
167   REAL, DIMENSION(:,:,:), ALLOCATABLE :: output_tab3d
[156]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
[1425]171   REAL,DIMENSION(:),ALLOCATABLE,SAVE :: dtrad
[156]172#endif
[11]173!-------------------------------------------
174! OUTPUT VARIABLES
175!-------------------------------------------
176
177
178!-------------------------------------------   
179! LOCAL  VARIABLES
180!-------------------------------------------
[170]181   INTEGER ::    i,j,k,its,ite,jts,jte,ij,kk
[73]182   INTEGER ::    subs,iii
[11]183
[170]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
[11]192   ! *** for LMD physics
193   ! ------> inputs:
[1038]194   INTEGER :: ngrid,nlayer,nq,nsoil
[11]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
[28]203#ifdef NEWPHYS
[678]204   REAL :: wstar_val,fluxrad_val
[28]205   REAL,DIMENSION(:),ALLOCATABLE :: isoil_val, dsoil_val
[315]206   REAL :: z0_val
[28]207#endif
[11]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 ?
[678]213   REAL,DIMENSION(:),ALLOCATABLE :: wtsurf,wco2ice,wemis
[34]214   REAL,DIMENSION(:,:),ALLOCATABLE :: wq2,wqsurf,wtsoil
[28]215#ifdef NEWPHYS
[678]216   REAL,DIMENSION(:),ALLOCATABLE :: wwstar,wfluxrad
[315]217   REAL,DIMENSION(:),ALLOCATABLE :: wz0tab
[28]218   REAL,DIMENSION(:,:),ALLOCATABLE :: wisoil,wdsoil
[55]219   CHARACTER*20,DIMENSION(:),ALLOCATABLE :: wtnom
[28]220#endif
[11]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, &
[55]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  :: &
[11]239     water_vapor_prof, water_ice_prof
[55]240#endif
[11]241
[55]242
[11]243   ! Additional control variables
244   INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe
[226]245   REAL :: elaps, ptimestep
[250]246   INTEGER :: wday_ini, test, test2
247   REAL :: wappel_phys
[11]248   LOGICAL :: flag_LES
[667]249   LOGICAL, SAVE :: flag_first_restart
[11]250!**************************************************
251! IMPORTANT: pour travailler avec grid nesting
[70]252! --- deux comportements distincts du save
253! ... ex1: ferme planeto avec PGF+MPI: standard
254! ... ex2: iDataPlex avec IFORT+MPI: SPECIAL_NEST_SAVE
[11]255!**************************************************
[70]256#ifdef SPECIAL_NEST_SAVE
257      !  une dimension supplementaire liee au nest
[67]258      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
[11]259             dp_save
[67]260      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
[11]261             du_save, dv_save, dt_save
[67]262      REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
[11]263             dq_save     
[688]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
[70]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     
[688]289#ifndef NORESTART
[674]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
[678]301#ifdef NEWPHYS
[674]302      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
303             save_wstar_restart
[678]304      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
305             save_fluxrad_restart
306#endif
[688]307#endif
308#endif
[674]309
[11]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
[34]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
[11]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!!
[94]361IF (flag_LES .eqv. .false.) THEN
[34]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
[11]366ips=its
367ipe=ite
368jps=jts
369jpe=jte
[94]370IF (flag_LES .eqv. .false.) THEN
[34]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
[11]376kps=kts         !! start at surface
[94]377IF (flag_LES .eqv. .false.) THEN
[34]378 kpe=kte-sponge_top
379ELSE
380 PRINT *, '*** IDEALIZED SIMULATION: LES *** kpe=kte'
381 kpe=kte !-sponge_top
382ENDIF
[11]383
384!----------------------------!
385! DIMENSIONS FOR THE PHYSICS !
386!----------------------------!
387wday_ini = JULDAY - 1      !! GCM convention
[250]388wappel_phys = RADT
389ptimestep = dt*wappel_phys            ! physical timestep (s)
[1212]390ngrid=(ipe-ips+1)*(jpe-jps+1)         ! size of the horizontal grid
[11]391nlayer = kpe-kps+1                    ! number of vertical layers: nlayermx
392nsoil = NUM_SOIL_LAYERS               ! number of soil layers: nsoilmx
[1038]393if (num_3d_s > 1) then                ! number of advected fields
[11]394        nq = num_3d_s-1               
395else
396        nq = 1
397endif
398! **** needed but hardcoded
399lastcall = .false.
400! **** needed but hardcoded
401
[155]402          !PRINT *, ips, ipe, jps, jpe
403          !PRINT *, ngrid
[11]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
[667]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
[11]416firstcall=.true.  !! for continuity with GCM, physics are always called at the first WRF timestep
[156]417  !firstcall=.false. !! just in case you'd want to get rid of the physics
[11]418test=0
[70]419#ifdef SPECIAL_NEST_SAVE
[67]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
[667]424   ALLOCATE(dp_save(ngrid,max_dom))
[67]425   ALLOCATE(du_save(ngrid,nlayer,max_dom))
[667]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))
[67]429   dp_save(:,:)=0.    !! initialize these arrays ...
430   du_save(:,:,:)=0.
431   dv_save(:,:,:)=0.
432   dt_save(:,:,:)=0.
433   dq_save(:,:,:,:)=0.
[688]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
[67]453ENDIF
[667]454IF (id .lt. max_dom) THEN
455   flag_first_restart=.true.
456ELSE
457   flag_first_restart=.false.
458ENDIF
[70]459#else
[700]460IF( .NOT. ALLOCATED( dp_save  ) ) THEN
[70]461ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
462ALLOCATE(du_save(ngrid,nlayer))
[667]463ALLOCATE(dv_save(ngrid,nlayer))
464ALLOCATE(dt_save(ngrid,nlayer))
465ALLOCATE(dq_save(ngrid,nlayer,nq))
[700]466ENDIF
[70]467dp_save(:)=0.    !! initialize these arrays ...
468du_save(:,:)=0.
469dv_save(:,:)=0.
[667]470dt_save(:,:)=0.
[70]471dq_save(:,:,:)=0.
[667]472flag_first_restart=.false.
[688]473#ifndef NORESTART
[674]474! Restart save arrays
[700]475IF( .NOT. ALLOCATED( save_tsoil_restart  ) ) THEN
[674]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))
[700]481ENDIF
[674]482save_tsoil_restart(:,:)=0.
483save_co2ice_restart(:)=0.
484save_q2_restart(:,:)=0.
485save_qsurf_restart(:,:)=0.
[678]486save_tsurf_restart(:)=0.
487#ifdef NEWPHYS
[700]488IF( .NOT. ALLOCATED( save_wstar_restart  ) ) THEN
[678]489ALLOCATE(save_wstar_restart(ngrid))
490ALLOCATE(save_fluxrad_restart(ngrid))
[700]491ENDIF
[674]492save_wstar_restart(:)=0.
[678]493save_fluxrad_restart(:)=0.
494#endif
[688]495#endif
496#endif
[674]497
[11]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
[667]500    print *, 'DOMAIN: ', ids, ide, jds, jde
[34]501    print *, 'MEMORY: ', ims, ime, jms, jme
[11]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!-------------------------------------------------!
[250]509IF (wappel_phys .gt. 0.) THEN
[11]510firstcall = .false.
[250]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)
[11]514ELSE
515firstcall = .false.
516test = 9999   
517ENDIF
518ENDIF
519
[341]520!!!!! for 'subgrid' temporal diagnostics
521!test2 = MODULO(elaps,history_interval*100.)
[11]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'
[34]542  PRINT *,'** Mars ** BEWARE: input_coord must be here'
[11]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)
[34]567!PRINT *, t
[11]568!        PRINT *,'************************************************'
569!        STOP
570!ENDIF   
571
[72]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)
[11]583
584!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
585!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[341]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.)
[11]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
[34]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)
[11]642
643
[34]644
[11]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
[70]662#ifdef SPECIAL_NEST_SAVE
[67]663pdpsrf(:)=dp_save(:,id)
664pdu(:,:)=du_save(:,:,id)
665pdv(:,:)=dv_save(:,:,id)
666pdt(:,:)=dt_save(:,:,id)
667pdq(:,:,:)=dq_save(:,:,:,id)
[70]668#else
669pdpsrf(:)=dp_save(:)
670pdu(:,:)=du_save(:,:)
671pdv(:,:)=dv_save(:,:)
672pdt(:,:)=dt_save(:,:)
673pdq(:,:,:)=dq_save(:,:,:)
674#endif
[11]675!!!
676!!! BIG LOOP : 2. calculate physical tendencies
677!!!
678ELSE
679!----------!
680! ALLOCATE !
681!----------!
682! inputs ...
[315]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))          !!!!!
[156]702#ifndef NOPHYS
[11]703ALLOCATE(output_tab2d(ngrid,n2d))
704ALLOCATE(output_tab3d(ngrid,nlayer,n3d))
[156]705#endif
[315]706ALLOCATE(wco2ice(ngrid))         !!!!!
707ALLOCATE(wemis(ngrid))           !!!!!
[11]708ALLOCATE(q2_val(nlayer+1))
709ALLOCATE(qsurf_val(nq))
710ALLOCATE(tsoil_val(nsoil))
[28]711#ifdef NEWPHYS
712ALLOCATE(isoil_val(nsoil))
713ALLOCATE(dsoil_val(nsoil))
714#endif
[315]715ALLOCATE(wq2(ngrid,nlayer+1))    !!!!!
716ALLOCATE(wqsurf(ngrid,nq))       !!!!!
717ALLOCATE(wtsoil(ngrid,nsoil))    !!!!!
[678]718#ifdef NEWPHYS
[674]719ALLOCATE(wfluxrad(ngrid))
720ALLOCATE(wwstar(ngrid))
[315]721ALLOCATE(wisoil(ngrid,nsoil))    !!!!!
722ALLOCATE(wdsoil(ngrid,nsoil))    !!!!!
[28]723#endif
[315]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))    !!!!!
[11]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))
[55]744#ifdef NEWPHYS
745ALLOCATE(q_prof(nlayer,nq))
[315]746ALLOCATE(wtnom(nq))              !!!!!
[55]747#else
[11]748ALLOCATE(water_vapor_prof(nlayer))
749ALLOCATE(water_ice_prof(nlayer))
[55]750#endif
[11]751
[315]752
[11]753!!!!!!!!!!!!!!!!!!!!!!!!!!!!
754!!!!!!!!!!!!!!!!!!!!!!!!!!!!
755!! PREPARE PHYSICS INPUTS !! 
756!!!!!!!!!!!!!!!!!!!!!!!!!!!!
757!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[76]758
[481]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
[76]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
[126]770!package   radioac      mars==20                     -              scalar:qtrac1
[170]771!package   radioac2     mars==21                     -              scalar:upward,downward
[324]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
[76]773!!!!!!!!!!!!!!!!!!! FOR REFERENCE
774
[55]775#ifdef NEWPHYS
[76]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 !!!!
[55]779SELECT CASE (MARS_MODE)
[802]780    CASE(0,10) 
[55]781      wtnom(nq) = 'co2'
[73]782    CASE(1)
[324]783      wtnom(1)  = 'h2o_vap'
784      wtnom(2)  = 'h2o_ice'
[73]785    CASE(2)
[76]786      wtnom(1)  = 'dust01'     
787    CASE(3)
788      wtnom(1)  = 'dust_mass'
[324]789      wtnom(2)  = 'dust_number'
[76]790    CASE(11)
791      wtnom(1)  = 'h2o_vap'
[324]792      wtnom(2)  = 'h2o_ice'
[76]793      wtnom(3)  = 'dust_mass'
794      wtnom(4)  = 'dust_number'
[481]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'
[126]802    CASE(20)
803      wtnom(1) = 'qtrac1'
[170]804    CASE(21)
805      wtnom(1) = 'upward'
[324]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'
[55]826END SELECT
827#endif
[11]828
[76]829!!*******************************************!!
830!!*******************************************!!
831
[11]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!--------------------------------!
[1388]863#ifndef NOPHYS
[55]864#ifdef NEWPHYS
[683]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
[802]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
[324]874  IF ((MARS_MODE .EQ. 20) .OR. (MARS_MODE .EQ. 21)) THEN
[674]875   IF (firstcall .EQV. .true. .and. (.not. restart)) THEN
[126]876      q_prof(:,:) = 0.95
[170]877   ENDIF
[324]878  ENDIF
[55]879#else
[11]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.
[76]890    CASE(3:)
891    print *, 'OOOPS... not ready yet.'
892    STOP
[11]893END SELECT
[55]894#endif
[1388]895#endif
[11]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!---------------------------------------------!
[34]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
[11]922
923!-----------------------------------------!
924! Gravity wave parametrization            !
925! NB: usually 0 in mesoscale applications !
926!-----------------------------------------!
[34]927IF (JULYR .ne. 9999) THEN
[1236]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)
[34]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
[11]941
942!---------------------------------!
943! Ground albedo & Thermal Inertia !
944!---------------------------------!
[34]945IF (JULYR .ne. 9999) THEN
946 IF (CST_AL == 0) THEN
[1236]947 albedodat_val=M_ALBEDO(i,j)
[34]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
[1236]953 inertiedat_val=M_TI(i,j)
[34]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
[11]963ENDIF
964
[315]965#ifdef NEWPHYS
966!----------------------------!
967! Variable surface roughness !
968!----------------------------!
969IF (JULYR .ne. 9999) THEN
970 IF (CST_Z0 == 0) THEN
[1236]971   z0_val = M_Z0(i,j)
[315]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
[324]982   IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
[315]983   z0_val = 0.01
984ENDIF
[1128]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
[315]992#endif
993
[11]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!-----------------------------------------------!
[1236]1003IF (M_TSURF(i,j) .gt. 0.) THEN
1004  tsurf_val=M_TSURF(i,j)
[674]1005ELSE
[1236]1006  tsurf_val=TSK(i,j) ! retro-compatibility
[674]1007ENDIF
[1236]1008emis_val=M_EMISS(i,j)
1009co2ice_val=M_CO2ICE(i,j)
[11]1010
1011!------------------------!
1012! Deep soil temperatures !
1013!------------------------!
[34]1014IF (JULYR .ne. 9999) THEN
[1236]1015  IF (M_TSOIL(i,1,j) .gt. 0.) THEN
1016   tsoil_val(:)=M_TSOIL(i,:,j)
[34]1017  ELSE
1018   tsoil_val = tsoil_val*0. + tsurf_val
1019  ENDIF
[54]1020#ifdef NEWPHYS
[1236]1021  isoil_val(:)=M_ISOIL(i,:,j)
1022  dsoil_val(:)=M_DSOIL(i,:,j)
[54]1023#endif
[11]1024ELSE
[34]1025  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION tsoil is set to tsurf'
1026  do k=1,nsoil
[674]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.
[1236]1031     !tsoil was saved in physiq.F under the name M_TSOIL in the restart file
[674]1032     !(see Registry)
[1236]1033     tsoil_val(k)=M_TSOIL(i,k,j)
[674]1034   ENDIF
1035
[54]1036#ifdef NEWPHYS
[574]1037   IF ( nsoil .lt. 18 ) THEN
[577]1038       PRINT *,'** Mars ** WRONG NUMBER OF SOIL LAYERS. SHOULD BE 18 AND IT IS ',nsoil
[574]1039       STOP
1040   ENDIF
[54]1041   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION isoil and dsoil standard'
1042   isoil_val(k) = inertiedat_val
[179]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
[54]1045   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** ISOIL DSOIL are ',isoil_val(k), dsoil_val(k)
1046#endif
[34]1047  enddo
[11]1048ENDIF
[34]1049
[11]1050!-------------------!
[72]1051! Tracer at surface !
1052!-------------------!
[1388]1053#ifndef NOPHYS
[72]1054SELECT CASE (MARS_MODE)
[76]1055    CASE(0) 
[72]1056    qsurf_val(:)=0.
[76]1057    CASE(1)
1058    qsurf_val(1)=0.
[1236]1059    qsurf_val(2)=M_H2OICE(i,j)    !! logique avec wtnom(2) = 'h2o_ice' defini ci-dessus
[76]1060                                   !! ----- retrocompatible ancienne physique
[1038]1061                                   !! ----- [H2O ice is last tracer in qsurf in LMD physics]
[76]1062    CASE(2) 
[77]1063    qsurf_val(1)=0.                !! not coupled with lifting for the moment [non remobilise]
[76]1064#ifdef NEWPHYS
1065    CASE(3)
[766]1066    qsurf_val(1)=q_prof(1,1)                !!! temporaire, a definir       
1067    qsurf_val(2)=q_prof(1,2)     
[485]1068    CASE(11)
[73]1069    qsurf_val(1)=0.
[1236]1070    qsurf_val(2)=M_H2OICE(i,j)    !! logique avec wtnom(2) = 'h2o_ice' defini ci-dessus
[77]1071    qsurf_val(3)=0.                !! not coupled with lifting for the moment [non remobilise]
[485]1072    qsurf_val(4)=0.
1073    CASE(12)   
1074    qsurf_val(1)=0.
[1236]1075    qsurf_val(2)=M_H2OICE(i,j)    !! logique avec wtnom(2) = 'h2o_ice' defini ci-dessus
[485]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.
[126]1080    CASE(20)
1081    qsurf_val(:)=0.
[170]1082    CASE(21)
1083    qsurf_val(:)=0.
[73]1084#else
[76]1085    CASE(3:)
1086    print *, 'OOOPS... not ready yet.'
1087    STOP
[73]1088#endif
[76]1089
[72]1090END SELECT
[1388]1091#endif
[72]1092
1093!-------------------!
[11]1094! Slope inclination !
1095!-------------------!
[34]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
[11]1103
1104!-------------------------------------------!
1105! Slope orientation; 0 is north, 90 is east !
1106!-------------------------------------------!
[34]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
[11]1119
1120!-----------------!
1121! Optional output !
1122!-----------------!
1123IF ( (i == ips) .AND. (j == jps) ) THEN
[34]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
[72]1140PRINT *,'qsurf ',qsurf_val
[28]1141#ifdef NEWPHYS
[34]1142PRINT *,'isoil ',isoil_val
1143PRINT *,'dsoil ',dsoil_val
[55]1144PRINT *,'q_prof ',q_prof
[315]1145PRINT *,'z0 ',z0_val
[28]1146#endif
[11]1147ENDIF
1148
1149!-------------------------!
1150!-------------------------!
1151!      PROVISOIRE         !
1152!-------------------------!
1153!-------------------------!
[674]1154IF (.not. restart) THEN
[700]1155   q2_val(:)=1.E-6      !PBL wind variance
[678]1156#ifdef NEWPHYS
[674]1157   fluxrad_val=0.
1158   wstar_val=0.
[678]1159#endif
[674]1160ELSE
[1236]1161   q2_val(:)=M_Q2(i,:,j)
[678]1162#ifdef NEWPHYS
[1236]1163   fluxrad_val=M_FLUXRAD(i,j)
1164   wstar_val=M_WSTAR(i,j)
[678]1165#endif
[674]1166ENDIF
[11]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(:)
[28]1188#ifdef NEWPHYS
[678]1189wfluxrad(subs) = fluxrad_val
1190wwstar(subs) = wstar_val
[28]1191wisoil(subs,:) = isoil_val(:)
1192wdsoil(subs,:) = dsoil_val(:)
[315]1193wz0tab(subs) = z0_val
[28]1194#endif
[11]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(:)
[28]1215pt(subs,:) = t_prof(:)
[11]1216pu(subs,:) = u_prof(:)
1217pv(subs,:) = v_prof(:)
1218pw(subs,:) = 0   !! NB: not used in the physics, only diagnostic...
[34]1219!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[11]1220!! for IDEALIZED CASES ONLY
1221IF (JULYR .eq. 9999) pplev(subs,nlayer+1)=0.  !! pplev(subs,nlayer+1)=ptop >> NO !
[34]1222!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[11]1223
[34]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
[11]1228!---------!
[55]1229! Tracers ! 
[11]1230!---------!
[1388]1231#ifndef NOPHYS
[55]1232#ifdef NEWPHYS
1233pq(subs,:,:) = q_prof(:,:)  !! traceurs generiques, seuls noms sont specifiques
1234#else
[11]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)
[76]1243    CASE(3:)
1244    print *, 'OOOPS... not ready yet.'
1245    STOP   
[11]1246END SELECT
[55]1247#endif
[1388]1248#endif
[11]1249ENDDO
1250ENDDO
1251
1252!!*****************!!
1253!! CLEAN THE PLACE !!
1254!!*****************!!
1255DEALLOCATE(q2_val)
1256DEALLOCATE(qsurf_val)
1257DEALLOCATE(tsoil_val)
[28]1258#ifdef NEWPHYS
1259DEALLOCATE(isoil_val)
1260DEALLOCATE(dsoil_val)
1261#endif
[11]1262DEALLOCATE(dz8w_prof)
1263DEALLOCATE(z_prof)
1264DEALLOCATE(p8w_prof)
1265DEALLOCATE(p_prof)
1266DEALLOCATE(t_prof)
[65]1267DEALLOCATE(t8w_prof)
[11]1268DEALLOCATE(u_prof)
1269DEALLOCATE(v_prof)
[55]1270#ifdef NEWPHYS
1271DEALLOCATE(q_prof)
1272#else
[11]1273DEALLOCATE(water_vapor_prof)
1274DEALLOCATE(water_ice_prof)
[55]1275#endif
[11]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
[156]1293#ifndef NOPHYS
[11]1294print *, '** Mars ** LMD INITIALIZATION'
[28]1295#include "../call_meso_inifis.inc"
1296!!! le # est important pour newphys
[156]1297#endif
[11]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)
[315]1311#ifdef NEWPHYS
1312DEALLOCATE(wz0tab)
1313#endif
[11]1314ENDIF
1315
1316!!********!!
1317!! PHYSIQ !!
1318!!********!!
[250]1319call_physics : IF (wappel_phys .ne. 0.) THEN
[11]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.
[156]1335#ifndef NOPHYS
1336print *, '** Mars ** CALL TO LMD PHYSICS'
[28]1337#include "../call_meso_physiq.inc"
1338!!! le # est important pour newphys
[156]1339#endif
[11]1340DEALLOCATE(pplev)
1341DEALLOCATE(pplay)
1342DEALLOCATE(pphi)
1343DEALLOCATE(pu)
1344DEALLOCATE(pv)
1345DEALLOCATE(pt)
1346DEALLOCATE(pw)
1347DEALLOCATE(pq)
[674]1348!DEALLOCATE(wtsurf)
1349!DEALLOCATE(wco2ice)
[11]1350DEALLOCATE(wemis)
[674]1351!DEALLOCATE(wq2)
1352!DEALLOCATE(wqsurf)
1353!DEALLOCATE(wtsoil)
1354!DEALLOCATE(wwstar)
1355!DEALLOCATE(wfluxrad)
[28]1356#ifdef NEWPHYS
1357DEALLOCATE(wisoil)
1358DEALLOCATE(wdsoil)
[55]1359DEALLOCATE(wtnom)
[28]1360#endif
[11]1361
1362
1363!-------------------------------!
1364! PHYSIQ OUTPUT IN THE WRF FILE !
1365!-------------------------------!
[156]1366#ifndef NOPHYS
[11]1367DO j = jps,jpe
1368DO i = ips,ipe
1369subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
[144]1370
[515]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
[144]1381
[11]1382#include "module_lmd_driver_output3.inc" 
1383       !  ^-- generated from Registry
1384ENDDO
1385ENDDO
1386DEALLOCATE(output_tab2d)
1387DEALLOCATE(output_tab3d)
[156]1388#endif
[11]1389
[118]1390!!!!!! PATCH SPECIAL STORM
1391#ifdef NEWPHYS
1392#ifdef storm
1393#include "../mars_lmd/storm.inc"
1394#endif
1395#endif
1396!!!!!! PATCH SPECIAL STORM
[117]1397
1398
[11]1399!---------------------------------------------------------------------------------!
1400! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
1401!---------------------------------------------------------------------------------!
[70]1402#ifdef SPECIAL_NEST_SAVE
[67]1403dp_save(:,id)=pdpsrf(:)
1404du_save(:,:,id)=pdu(:,:)
1405dv_save(:,:,id)=pdv(:,:)
1406dt_save(:,:,id)=pdt(:,:)
1407dq_save(:,:,:,id)=pdq(:,:,:)
[688]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
[70]1419#else
1420dp_save(:)=pdpsrf(:)
1421du_save(:,:)=pdu(:,:)
1422dv_save(:,:)=pdv(:,:)
1423dt_save(:,:)=pdt(:,:)
1424dq_save(:,:,:)=pdq(:,:,:)
[688]1425#ifndef NORESTART
[674]1426save_tsoil_restart(:,:)=wtsoil(:,:)
1427save_co2ice_restart(:)=wco2ice(:)
1428save_q2_restart(:,:)=wq2(:,:)
1429save_qsurf_restart(:,:)=wqsurf(:,:)
1430save_tsurf_restart(:)=wtsurf(:)
[688]1431#ifdef NEWPHYS
1432save_wstar_restart(:)=wwstar(:)
1433save_fluxrad_restart(:)=wfluxrad(:)
1434#endif
1435#endif
1436#endif
[674]1437DEALLOCATE(wtsoil)
1438DEALLOCATE(wco2ice)
1439DEALLOCATE(wq2)
1440DEALLOCATE(wqsurf)
[678]1441DEALLOCATE(wtsurf)
1442#ifdef NEWPHYS
1443DEALLOCATE(wfluxrad)
[674]1444DEALLOCATE(wwstar)
[678]1445#endif
[11]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!------------------------------------------------------------------!
[170]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'
[503]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
[170]1498ENDIF
1499
[1388]1500#ifdef NOPHYS
1501!!!!!!!!!!!!!!!!!!!
1502!!!!!!!!!!!!!!!!!!!
1503IF (firstcall .EQV. .true.) THEN
[1427]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)
[1388]1511ENDIF
1512DO i= 1,ngrid
[1427]1513  pdt(i,kps:kpe) = dtrad(kps:kpe)
[1388]1514ENDDO
1515!!!!!!!!!!!!!!!!!!!
1516!!!!!!!!!!!!!!!!!!!
1517#endif
1518
[11]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!---------------------------!
[117]1543PSFC(i,j)=PSFC(i,j)+pdpsrf(subs)*dt    !!! here dt is needed
[11]1544
[674]1545!------------------------------------!
1546! Save key variables for restart ! 
1547!------------------------------------!
[688]1548#ifndef NORESTART
1549#ifdef SPECIAL_NEST_SAVE
[1236]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)
[688]1553SELECT CASE (MARS_MODE)
1554   CASE (1,11,12)
[1236]1555     M_H2OICE(i,j)=save_qsurf_restart(subs,2,id)  !! see above Tracer at surface
[688]1556END SELECT
[1236]1557M_TSURF(i,j)=save_tsurf_restart(subs,id)
[688]1558#ifdef NEWPHYS
[1236]1559M_WSTAR(i,j)=save_wstar_restart(subs,id)
1560M_FLUXRAD(i,j)=save_fluxrad_restart(subs,id)
[688]1561#endif
1562#else
[1236]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,:)
[674]1566SELECT CASE (MARS_MODE)
1567   CASE (1,11,12)
[1236]1568     M_H2OICE(i,j)=save_qsurf_restart(subs,2)  !! see above Tracer at surface
[674]1569END SELECT
[1236]1570M_TSURF(i,j)=save_tsurf_restart(subs)
[678]1571#ifdef NEWPHYS
[1236]1572M_WSTAR(i,j)=save_wstar_restart(subs)
1573M_FLUXRAD(i,j)=save_fluxrad_restart(subs)
[678]1574#endif
[688]1575#endif
1576#endif
[674]1577
[11]1578!---------!
[55]1579! Tracers ! 
[11]1580!---------!
[1388]1581#ifndef NOPHYS
[55]1582#ifdef NEWPHYS
1583SCALAR(i,kps:kpe,j,1)=0.
[170]1584
1585SELECT CASE (MARS_MODE)
[683]1586CASE(0)
1587        SCALAR(i,kps:kpe,j,:)=0.
[170]1588CASE(20)
1589   !! Mars mode 20 : add a passive tracer with radioactive-like decay
[129]1590   IF ( (i == ips) .AND. (j == jps) )   print *, 'RADIOACTIVE-LIKE TRACER WITH SOURCE AT SURFACE LAYER.'
[126]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
[170]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
[126]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
[170]1607END SELECT
[55]1608#else
[11]1609SELECT CASE (MARS_MODE)
1610CASE(0)
1611        SCALAR(i,kps:kpe,j,:)=0.
1612CASE(1)
1613        !!! Water vapor
[117]1614        SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*dt  !!! here dt is needed
[11]1615        !!! Water ice
[117]1616        SCALAR(i,kps:kpe,j,3)=SCALAR(i,kps:kpe,j,3)+pdq(subs,kps:kpe,nq-1)*dt  !!! here dt is needed
[11]1617CASE(2)
1618        !!! Dust
[117]1619        SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*dt  !!! here dt is needed
[76]1620CASE(3:)
1621        print *, 'OOOPS... not ready yet.'
1622        STOP   
[11]1623END SELECT
[55]1624#endif
[1388]1625#endif
[11]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
[34]1655!!---------!
1656!! display !
1657!!---------!
[155]1658!PRINT *, '** Mars ** Results from LMD physics'
1659!PRINT *, 'u non-zero tendencies'
[67]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.)
[72]1664!PRINT *, RUBLTEN(10,1,10), RUBLTEN(10,15,10)
[155]1665!PRINT *, 'v non-zero tendencies'
[67]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.)
[72]1670!PRINT *, RVBLTEN(10,1,10), RVBLTEN(10,15,10)
[11]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
[34]1697!c      xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
[11]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.