source: trunk/WRF.COMMON/WRFV3/modif_mars/module_lmd_driver.F.bak @ 2757

Last change on this file since 2757 was 44, checked in by aslmd, 14 years ago

LMD_MM_MARS: correction mineure (fichiers manquants sur le serveur)
File size: 31.1 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!*******************************************************************************
12MODULE module_lmd_driver
13CONTAINS
14    SUBROUTINE lmd_driver(id,max_dom,DT,ITIMESTEP,XLAT,XLONG, &
15        IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME,KMS,KME, &
16        i_start,i_end,j_start,j_end,kts,kte,num_tiles, &
17        DX,DY, &
18        MSFT,MSFU,MSFV, &
19        GMT,JULYR,JULDAY, &
20        P8W,DZ8W,T8W,Z,HT, &
21        U,V,TH,T,P,EXNER,RHO, &
22        PTOP, &
23        RADT,CUDT, &
24        TSK,PSFC, &
25        RTHBLTEN,RUBLTEN,RVBLTEN, &
26        num_3d_s,SCALAR, &
27        MARS_MODE, &
28        MARS_ALB,MARS_TI,MARS_CICE,MARS_EMISS, &
29        MARS_TSOIL, &
30        MARS_GW, &
31        NUM_SOIL_LAYERS, &
32        CST_AL, &
33        CST_TI, &
34        isfflx, &
35#include "../modif_mars/module_lmd_driver_output1.inc"
36        SLPX,SLPY)
37! NB: module_lmd_driver_output1.inc : output arguments generated from Registry
38
39!        IPS,IPE,JPS,JPE,KPS,KPE, &
40
41
42
43!==================================================================
44! USES
45!==================================================================
46   USE module_model_constants
47   USE module_wrf_error
48! add new modules here, if needed ...
49
50!==================================================================
51  IMPLICIT NONE
52!==================================================================
53
54!==================================================================
55! COMMON
56!==================================================================
57
58!! the only common needed is the one defining the physical grid
59!! -- path is hardcoded, but the structure is not subject to change
60!! -- please put # if needed by the pre-compilation process     
61!
62!
63include "../modif_mars/dimphys.h"
64!
65!--to be commented because there are tests in the physics ?
66!--TODO: get rid of the ...mx first in this routine and .inc
67
68!
69! INCLUDE AUTOMATIQUEMENT GENERE A PARTIR DU REGISTRY
70!
71include "../modif_mars/wrf_output_2d.h"
72include "../modif_mars/wrf_output_3d.h"
73
74!==================================================================
75! VARIABLES
76!==================================================================
77
78! WRF Dimensions
79INTEGER, INTENT(IN   )    ::   &
80      ids,ide,jds,jde,kds,kde, &
81      ims,ime,jms,jme,kms,kme, &
82!      ips,ipe,jps,jpe,kps,kpe, &
83            kts,kte,num_tiles, &
84              NUM_SOIL_LAYERS
85INTEGER, DIMENSION(num_tiles), INTENT(IN) ::  &
86      i_start,i_end,j_start,j_end   
87! Scalars
88INTEGER, INTENT(IN  ) :: JULDAY, itimestep, julyr,id,max_dom,MARS_MODE,isfflx
89REAL, INTENT(IN  ) :: GMT,dt,dx,dy,RADT,CUDT
90REAL, INTENT(IN  ) :: CST_AL, CST_TI
91REAL, INTENT(IN  ) :: PTOP
92! 2D arrays         
93REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  :: &
94     MSFT,MSFU,MSFV, &
95     XLAT,XLONG,HT,  &
96     MARS_ALB,MARS_TI,MARS_EMISS,MARS_CICE, &
97     SLPX,SLPY
98! 3D arrays
99REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
100     dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th
101REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(IN   )  :: &
102     MARS_TSOIL
103REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: &
104     MARS_GW
105! 4D arrays
106INTEGER, INTENT(IN ) :: num_3d_s
107REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_s ), INTENT(INOUT ) :: &
108     scalar
109
110!-------------------------------------------
111! OUTPUT VARIABLES
112!-------------------------------------------
113!
114! Generated from Registry
115!
116! default definitions :
117! 2D : TSK, PSFC
118! 3D : RTHBLTEN,RUBLTEN,RVBLTEN
119#include "../modif_mars/module_lmd_driver_output2.inc"
120   REAL, DIMENSION(:,:), ALLOCATABLE :: output_tab2d
121   REAL, DIMENSION(:,:,:), ALLOCATABLE :: output_tab3d
122!-------------------------------------------
123! OUTPUT VARIABLES
124!-------------------------------------------
125
126
127!-------------------------------------------   
128! LOCAL  VARIABLES
129!-------------------------------------------
130   INTEGER ::    i,j,k,ij 
131   INTEGER ::    its,ite,jts,jte
132   INTEGER ::    subs
133
134   ! *** for LMD physics
135   ! ------> inputs:
136   INTEGER :: ngrid,nlayer,nq,nsoil,nqmx
137   REAL :: pday,ptime,MY 
138   REAL :: aire_val,lat_val,lon_val
139   REAL :: phisfi_val,albedodat_val,inertiedat_val
140   REAL :: tsurf_val,co2ice_val,emis_val
141   REAL :: zmea_val,zstd_val,zsig_val,zgam_val,zthe_val
142   REAL :: theta_val, psi_val
143   LOGICAL :: firstcall,lastcall,tracerdyn
144   REAL,DIMENSION(:),ALLOCATABLE :: q2_val, qsurf_val, tsoil_val
145   REAL,DIMENSION(:),ALLOCATABLE :: aire_vec,lat_vec,lon_vec
146   REAL,DIMENSION(:),ALLOCATABLE :: walbedodat,winertiedat,wphisfi
147   REAL,DIMENSION(:),ALLOCATABLE :: wzmea,wzstd,wzsig,wzgam,wzthe
148   REAL,DIMENSION(:),ALLOCATABLE :: wtheta, wpsi
149   ! v--- can they be modified ?
150   REAL,DIMENSION(:),ALLOCATABLE :: wtsurf,wco2ice,wemis
151   REAL,DIMENSION(:,:),ALLOCATABLE :: wq2,wqsurf,wtsoil
152   ! ----------
153   REAL,DIMENSION(:,:),ALLOCATABLE :: pplev,pplay,pphi,pu,pv,pt,pw
154   REAL,DIMENSION(:,:,:),ALLOCATABLE :: pq 
155
156   ! <------ outputs:
157   !     physical tendencies
158   REAL,DIMENSION(:),ALLOCATABLE :: pdpsrf
159   REAL,DIMENSION(:,:),ALLOCATABLE :: pdu,pdv,pdt
160   REAL,DIMENSION(:,:,:),ALLOCATABLE :: pdq
161   ! ... intermediate arrays
162   REAL, DIMENSION(:), ALLOCATABLE  :: &
163     dz8w_prof,p8w_prof,p_prof,t_prof,t8w_prof, &
164     u_prof,v_prof,z_prof, &
165     water_vapor_prof, water_ice_prof
166!!   pi_prof, rho_prof, th_prof, &
167
168   ! Additional control variables
169   INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe
170   REAL :: elaps, ptimestep, wecri_phys_sec
171   INTEGER :: wappel_phys, wecri_phys, wday_ini, test
172   LOGICAL :: flag_LES
173
174!**************************************************
175! IMPORTANT: pour travailler avec grid nesting
176!       mieux vaut ne pas utiliser de SAVE ??
177!**************************************************
178      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
179             dp_save
180      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
181             du_save, dv_save, dt_save
182      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
183             dq_save     
184
185!!!IDEALIZED IDEALIZED
186      REAL :: lat_input, lon_input, ls_input, lct_input
187!!!IDEALIZED IDEALIZED
188
189
190
191!==================================================================
192! CODE
193!==================================================================
194
195!! SPECIAL LES
196IF (isfflx .ne. 0) THEN
197     flag_LES = .true.
198ELSE
199     flag_LES = .false.
200ENDIF
201!! SPECIAL LES
202
203print *,'** Mars ** DOMAIN',id
204
205!-------------------------!
206! TWEAK on WRF DIMENSIONS !
207!-------------------------!
208its = i_start(1)   ! define here one big tile
209ite = i_end(num_tiles)
210jts = j_start(1)
211jte = j_end(num_tiles)
212!!
213!relax=0
214!sponge_top=0               ! another value than 0 triggers instabilities 
215!IF (id .gt. 1) relax=2     ! fix to avoid noise in nesting simulations ; 1 >> too much noise ...
216ips=its
217ipe=ite
218jps=jts
219jpe=jte
220!IF (ips .eq. ids)   ips=its+relax !! IF tests necesary for parallel runs
221!IF (ipe .eq. ide-1) ipe=ite-relax
222!IF (jps .eq. jds)   jps=jts+relax
223!IF (jpe .eq. jde-1) jpe=jte-relax
224kps=kts         !! start at surface
225kpe=kte !-sponge_top
226
227!----------------------------!
228! DIMENSIONS FOR THE PHYSICS !
229!----------------------------!
230wday_ini = JULDAY - 1      !! GCM convention
231wappel_phys = int(RADT)
232wecri_phys = int(CUDT)
233ptimestep = dt*float(wappel_phys)     ! physical timestep (s)
234ngrid=(ipe-ips+1)*(jpe-jps+1)         ! size of the horizontal grid: ngridmx = wiim * wjjm
235nlayer = kpe-kps+1                    ! number of vertical layers: nlayermx
236nsoil = NUM_SOIL_LAYERS               ! number of soil layers: nsoilmx
237if (num_3d_s > 1) then                ! number of advected fields: nqmx
238        nq = num_3d_s-1               
239        nqmx = num_3d_s-1
240else
241        nq = 1
242        nqmx = 1
243endif
244! **** needed but hardcoded
245lastcall = .false.
246! **** needed but hardcoded
247
248
249PRINT *, ips, ipe, jps, jpe
250PRINT *, ngrid
251
252
253
254elaps = (float(itimestep)-1.)*dt  ! elapsed seconds of simulation
255!----------------------------------------------!
256! what is done at the first step of simulation !
257!----------------------------------------------!
258IF (elaps .eq. 0.) THEN
259firstcall=.true.  !! for continuity with GCM, physics are always called at the first WRF timestep
260test=0
261ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
262ALLOCATE(du_save(ngrid,nlayer))
263ALLOCATE(dv_save(ngrid,nlayer))
264ALLOCATE(dt_save(ngrid,nlayer))     
265ALLOCATE(dq_save(ngrid,nlayer,nq))
266dp_save(:)=0.    !! initialize these arrays ...
267du_save(:,:)=0.
268dv_save(:,:)=0.
269dt_save(:,:)=0.   
270dq_save(:,:,:)=0.
271!! put here some general information you'd like to print just once
272    print *, 'TILES: ', i_start,i_end, j_start, j_end  ! numbers for simple runs, arrays for parallel runs
273    print *, 'DOMAIN: ', ids, ide, jds, jde
274    print *, 'MEMORY: ', ims, ime, jms, jme
275    print *, 'ADVECTED TRACERS: ', num_3d_s-1
276    print *, 'PHYSICS IS CALLED EACH...',wappel_phys
277!! put here some general information you'd like to print just once
278ELSE
279!-------------------------------------------------!
280! what is done for the other steps of simulation  !
281!-------------------------------------------------!
282IF (wappel_phys .gt. 0) THEN
283firstcall = .false.
284test = MODULO(itimestep-1,wappel_phys) ! WRF time is integrated time (itimestep=1 at the end of first step)
285                                       ! -- same strategy as diagfi in the LMD GCM
286ELSE
287firstcall = .false.
288test = 9999   
289ENDIF
290ENDIF
291
292
293!!!******!!
294!!! TIME !!
295!!!******!!
296IF (JULYR .ne. 9999) THEN
297  !
298  ! specified
299  !
300  ptime = (GMT + elaps/3700.) !! universal time (0<ptime<1): ptime=0.5 at 12:00 UT
301  ptime = MODULO(ptime,24.)   !! the two arguments of MODULO must be of the same type
302  ptime = ptime / 24.
303  pday = (JULDAY - 1 + INT((3700*GMT+elaps)/88800))
304  pday = MODULO(int(pday),669)
305  MY = (JULYR-2000) + (88800.*(JULDAY - 1)+3700.*GMT+elaps)/59496000.
306  MY = INT(MY)
307ELSE
308  !
309  ! idealized
310  !
311  PRINT *,'** Mars ** IDEALIZED SIMULATION'
312  open(unit=14,file='input_coord',form='formatted',status='old')
313  rewind(14)
314  read(14,*) lon_input
315  read(14,*) lat_input
316  read(14,*) ls_input
317  read(14,*) lct_input
318  close(14)
319  ptime = lct_input - lon_input / 15. + elaps/3700.
320  ptime = MODULO(ptime,24.)
321  ptime = ptime / 24.
322  pday = floor(ls2sol(ls_input)) + INT((3700*(lct_input - lon_input / 15.) + elaps)/88800)
323  pday = MODULO(int(pday),669)
324  MY = 2024
325  wday_ini = floor(ls2sol(ls_input))
326ENDIF
327print *,'** Mars ** TIME IS', pday, ptime*24.
328
329
330!!****************!!
331!! CHECK DYNAMICS !!
332!!****************!!
333!IF ((MAXVAL(t) > 500).OR.(MINVAL(t,MASK = t > 0) <= 50)) THEN
334!        PRINT *,'******************   CRASH   *******************'
335!        PRINT *,'Irrealistic temperature...', MAXLOC(t), MINLOC(t)
336!PRINT *, t
337!        PRINT *,'************************************************'
338!        STOP
339!ENDIF   
340!
341!!! PB SAUF SI debug = 200 ?!!?
342!IF (float(itimestep) > 200.) THEN        ! to allow initialisation with zero-wind or constant
343!                                         ! and allow some outputs to locate the NaNs :)
344!IF (         (abs(MAXVAL(u)) == 0.) &
345!        .OR. (MINVAL(u) > MAXVAL(u)) &
346!        .OR. (abs(MAXVAL(v)) == 0.) &
347!        .OR. (MINVAL(v) > MAXVAL(v)) ) THEN
348!!IF (          (ANY(isNaN(u)) .EQV. .true.) &
349!!         .OR. (ANY(isNaN(v)) .EQV. .true.) & 
350!!         .OR. (ANY(isNaN(t)) .EQV. .true.) ) THEN
351!IF (         ANY(u*0. /= 0.) &
352!        .OR. ANY(v*0. /= 0.)        ) THEN
353!        PRINT *,'******************   CRASH   *******************'
354!        PRINT *,'************************************************'
355!        PRINT *,'NaN appeared in the simulation ...'
356!        PRINT *,'...this may be due to numerical or dynamical instability'
357!        PRINT *,'************************************************'
358!        PRINT *,'POSSIBLE SOLUTIONS:'
359!        PRINT *,'>> IF nonhydrostatic mode,'
360!        PRINT *,'   --> check that smdiv, emdiv and epssm are not 0.'
361!        PRINT *,'>> IF cfl is violated, '
362!        PRINT *,'   --> try to lower the dynamical timestep'
363!        PRINT *,'>> IF topographical gradients are high near specified bdy,'
364!        PRINT *,'   --> try to redefine the domain'
365!        PRINT *,'************************************************'
366!        STOP
367!ENDIF
368!ENDIF
369        !IF (          ANY(isNaN(u)) &
370        !         .OR. ANY(isNaN(v)) &
371        !         .OR. ANY(isNaN(t)) ) THEN
372        ! >>> ne marche qu'avec g95
373!print *, 'check dynamics'
374!print *, 'u', MAXVAL(u), MINVAL(u)
375!print *, 'v', MAXVAL(v), MINVAL(v)
376!print *, 't', MAXVAL(t), MINVAL(t, MASK = t > 0)
377
378
379
380!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
381!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
382!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
383!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
384!----------!
385! ALLOCATE !
386!----------!
387ALLOCATE(pdpsrf(ngrid))
388ALLOCATE(pdu(ngrid,nlayer))
389ALLOCATE(pdv(ngrid,nlayer))
390ALLOCATE(pdt(ngrid,nlayer))
391ALLOCATE(pdq(ngrid,nlayer,nq))
392!!!
393!!! BIG LOOP : 1. no call for physics, used saved values
394!!!
395IF (test.NE.0) THEN
396print *,'** Mars ** NO CALL FOR PHYSICS, go to next step...',test
397pdpsrf(:)=dp_save(:)
398pdu(:,:)=du_save(:,:)
399pdv(:,:)=dv_save(:,:)
400pdt(:,:)=dt_save(:,:)
401pdq(:,:,:)=dq_save(:,:,:)
402        !!!!!TEST TEST
403        !!!!!pour le nesting c est mieux de les mettre dans la physique, les SAVE
404!!!
405!!! BIG LOOP : 2. calculate physical tendencies
406!!!
407ELSE
408!----------!
409! ALLOCATE !
410!----------!
411! inputs ...
412ALLOCATE(aire_vec(ngrid))
413ALLOCATE(lon_vec(ngrid))
414ALLOCATE(lat_vec(ngrid))
415ALLOCATE(walbedodat(ngrid))
416ALLOCATE(winertiedat(ngrid))
417ALLOCATE(wphisfi(ngrid))
418ALLOCATE(wzmea(ngrid))
419ALLOCATE(wzstd(ngrid))
420ALLOCATE(wzsig(ngrid))
421ALLOCATE(wzgam(ngrid))
422ALLOCATE(wzthe(ngrid))
423ALLOCATE(wtheta(ngrid))
424ALLOCATE(wpsi(ngrid))
425ALLOCATE(wtsurf(ngrid))
426ALLOCATE(output_tab2d(ngrid,n2d))
427ALLOCATE(output_tab3d(ngrid,nlayer,n3d))
428ALLOCATE(wco2ice(ngrid))
429ALLOCATE(wemis(ngrid))
430ALLOCATE(q2_val(nlayer+1))
431ALLOCATE(qsurf_val(nq))
432ALLOCATE(tsoil_val(nsoil))
433ALLOCATE(wq2(ngrid,nlayer+1))
434ALLOCATE(wqsurf(ngrid,nq))
435ALLOCATE(wtsoil(ngrid,nsoil))
436ALLOCATE(pplev(ngrid,nlayer+1))
437ALLOCATE(pplay(ngrid,nlayer))
438ALLOCATE(pphi(ngrid,nlayer))
439ALLOCATE(pu(ngrid,nlayer))
440ALLOCATE(pv(ngrid,nlayer))
441ALLOCATE(pt(ngrid,nlayer))
442ALLOCATE(pw(ngrid,nlayer))
443ALLOCATE(pq(ngrid,nlayer,nq))
444! interm
445ALLOCATE(dz8w_prof(nlayer))
446ALLOCATE(p8w_prof(nlayer))
447ALLOCATE(p_prof(nlayer))
448ALLOCATE(t_prof(nlayer))
449ALLOCATE(t8w_prof(nlayer))
450ALLOCATE(u_prof(nlayer))
451ALLOCATE(v_prof(nlayer))
452ALLOCATE(z_prof(nlayer))
453!ALLOCATE(th_prof(nlayer))
454!ALLOCATE(rho_prof(nlayer))
455!ALLOCATE(pi_prof(nlayer))
456ALLOCATE(water_vapor_prof(nlayer))
457ALLOCATE(water_ice_prof(nlayer))
458
459
460!!!!!!!!!!!!!!!!!!!!!!!!!!!!
461!!!!!!!!!!!!!!!!!!!!!!!!!!!!
462!! PREPARE PHYSICS INPUTS !! 
463!!!!!!!!!!!!!!!!!!!!!!!!!!!!
464!!!!!!!!!!!!!!!!!!!!!!!!!!!!
465
466
467DO j = jps,jpe
468DO i = ips,ipe
469
470!!*******************************************!!
471!! FROM 3D WRF FIELDS TO 1D PHYSICS PROFILES !!
472!!*******************************************!!
473
474!-----------------------------------!
475! 1D subscript for physics "cursor" !
476!-----------------------------------!
477subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
478
479!--------------------------------------!
480! 1D columns : inputs for the physics  !
481! ... vert. coord., meteor. fields ... !
482!--------------------------------------!
483dz8w_prof(:) = dz8w(i,kps:kpe,j)   ! dz between full levels (m)   
484p8w_prof(:) = p8w(i,kps:kpe,j)     ! pressure full level (Pa) >> pplev
485p_prof(:) = p(i,kps:kpe,j)         ! pressure half level (Pa) >> pplay
486t_prof(:) = t(i,kps:kpe,j)         ! temperature half level (K) >> pt
487t8w_prof(:) = t8w(i,kps:kpe,j)     ! temperature full level (K)
488u_prof(:) = u(i,kps:kpe,j)         ! zonal wind (A-grid: unstaggered) half level (m/s) >> pu 
489v_prof(:) = v(i,kps:kpe,j)         ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv
490z_prof(:) = z(i,kps:kpe,j)         ! geopotential height half level (m) >> pphi/g
491!pi_prof(:) = exner(i,kps:kpe,j)    ! exner function (dimensionless) half level
492!rho_prof(:) = rho(i,kps:kpe,j)     ! density half level 
493!th_prof(:) = th(i,kps:kpe,j)       ! pot. temperature half level (K)
494
495!--------------------------------!
496! specific treatment for tracers !
497!--------------------------------!
498SELECT CASE (MARS_MODE)
499    CASE(0)  !! NO TRACERS (mars=0)
500    water_vapor_prof(:) = 0.
501    water_ice_prof(:) = 0.
502    CASE(1)  !! WATER CYCLE (mars=1)
503    water_vapor_prof(:) = scalar(i,kps:kpe,j,2)  !! H2O vapor is tracer 1 in the Registry for mars=1
504    water_ice_prof(:) = scalar(i,kps:kpe,j,3)    !! H2O ice is tracer 2 in the Registry for mars=1
505    CASE(2)  !! DUST CYCLE (mars=2)
506    water_vapor_prof(:) = 0.
507    water_ice_prof(:) = 0.
508END SELECT
509
510!!**********************************************************!!
511!! STATIC FIELDS FOR THE PHYSICS - NEEDED ONLY AT FIRSTCALL !!
512!!**********************************************************!!
513needed_ini_phys : IF (firstcall .EQV. .true.) THEN
514
515!----------------------------------------!
516! Surface of each part of the grid (m^2) !
517!----------------------------------------!
518!aire_val = dx*dy                           !! 1. idealized cases - computational grid
519aire_val = (dx/msft(i,j))*(dy/msft(i,j))    !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
520!aire_val=dx*dy/msfu(i,j)                   !! 3. special for Mercator GCM-like simulations
521
522!!---------------------------------------------!
523!! Mass-point latitude and longitude (radians) !
524!!---------------------------------------------!
525!lat_val = XLAT(i,j)*DEGRAD
526!lon_val = XLONG(i,j)*DEGRAD
527lat_val = lat_input*DEGRAD
528lon_val = lon_input*DEGRAD
529
530!!-----------------------------------------!
531!! Gravity wave parametrization            !
532!! NB: usually 0 in mesoscale applications !
533!!-----------------------------------------!
534!zmea_val=MARS_GW(i,1,j)
535!zstd_val=MARS_GW(i,2,j)
536!zsig_val=MARS_GW(i,3,j)
537!zgam_val=MARS_GW(i,4,j)
538!zthe_val=MARS_GW(i,5,j)
539zmea_val=0.
540zstd_val=0.
541zsig_val=0.
542zgam_val=0.
543zthe_val=0.
544
545!!---------------------------------!
546!! Ground albedo & Thermal Inertia !
547!!---------------------------------!
548!albedodat_val=MARS_ALB(i,j)
549!inertiedat_val=MARS_TI(i,j)
550albedodat_val=CST_AL
551inertiedat_val=CST_TI
552
553!---------------------------------------------!
554! Ground geopotential                         !
555! (=g*HT(i,j), only used in the microphysics) !
556!---------------------------------------------!
557phisfi_val=g*(z_prof(1)-dz8w_prof(1)/2.)   !! NB: z_prof are half levels, so z_prof(1) is not the surface
558
559!-----------------------------------------------!
560! Ground temperature, emissivity, CO2 ice cover !
561!-----------------------------------------------!
562tsurf_val=tsk(i,j)
563emis_val=MARS_EMISS(i,j)
564co2ice_val=MARS_CICE(i,j)
565
566!!------------------------!
567!! Deep soil temperatures !
568!!------------------------!
569!tsoil_val(:)=MARS_TSOIL(i,:,j)
570do k=1,nsoil
571 tsoil_val(k) = tsurf_val
572enddo
573
574!!-------------------!
575!! Slope inclination !
576!!-------------------!
577!theta_val=atan(sqrt( (1000.*SLPX(i,j))**2 + (1000.*SLPY(i,j))**2 ))
578!theta_val=theta_val/DEGRAD
579theta_val=0.
580
581!!-------------------------------------------!
582!! Slope orientation; 0 is north, 90 is east !
583!!-------------------------------------------!
584!psi_val=-90.*DEGRAD-atan(SLPY(i,j)/SLPX(i,j))
585!if (SLPX(i,j) .ge. 0.) then
586!   psi_val=psi_val-180.*DEGRAD
587!endif
588!psi_val=360.*DEGRAD+psi_val
589!psi_val=psi_val/DEGRAD
590!psi_val = MODULO(psi_val+180.,360.)
591psi_val=0.
592
593!
594! PRINT
595!
596IF ( (i == ips) .AND. (j == jps) ) THEN
597PRINT *,'lat/lon ', lat_val/DEGRAD, lon_val/DEGRAD
598PRINT *,'emiss ', emis_val
599PRINT *,'albedo ', albedodat_val
600PRINT *,'inertie ', inertiedat_val
601PRINT *,'phi ',phisfi_val
602PRINT *,'tsurf ',tsurf_val
603PRINT *,'aire ',aire_val
604PRINT *,'z_prof ',z_prof
605PRINT *,'dz8w_prof',dz8w_prof
606PRINT *,'p8w_prof ',p8w_prof
607PRINT *,'p_prof ',p_prof
608PRINT *,'t_prof ',t_prof
609PRINT *,'t8w_prof ',t8w_prof
610PRINT *,'u_prof ',u_prof
611PRINT *,'v_prof ',v_prof
612PRINT *,'tsoil ',tsoil_val
613ENDIF
614
615!-------------------------!
616!-------------------------!
617!      PROVISOIRE         !
618!-------------------------!
619!-------------------------!
620q2_val(:)=0      !PBL wind variance
621qsurf_val(:)=0   !Tracer on surface
622
623!-----------------!
624! Fill the arrays !
625!-----------------!
626aire_vec(subs) = aire_val   !! NB: total area in square km is SUM(aire_vec)/1.0E6
627lat_vec(subs) = lat_val
628lon_vec(subs) = lon_val
629wphisfi(subs) = phisfi_val
630walbedodat(subs) = albedodat_val
631winertiedat(subs) = inertiedat_val
632wzmea(subs) = zmea_val
633wzstd(subs) = zstd_val
634wzsig(subs) = zsig_val
635wzgam(subs) = zgam_val
636wzthe(subs) = zthe_val
637wtsurf(subs) = tsurf_val
638wco2ice(subs) = co2ice_val
639wemis(subs) = emis_val
640wq2(subs,:) = q2_val(:)
641wqsurf(subs,:) = qsurf_val(:)
642wtsoil(subs,:) = tsoil_val(:)
643wtheta(subs) = theta_val
644wpsi(subs) = psi_val
645
646ENDIF needed_ini_phys
647
648!!*****************************!!
649!! PREPARE "FOOD" FOR PHYSIQ.F !!
650!!*****************************!!
651
652!---------------------------------------------!
653! in LMD physics, geopotential must be        !
654! expressed with respect to the local surface !
655!---------------------------------------------!
656pphi(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) )
657
658!--------------------------------!
659! Dynamic fields for LMD physics !
660!--------------------------------!
661pplev(subs,1:nlayer) = p8w_prof(1:nlayer)  !! NB: last level: no data
662pplay(subs,:) = p_prof(:)
663pt(subs,:) = t_prof(:)
664pu(subs,:) = u_prof(:)
665pv(subs,:) = v_prof(:)
666pw(subs,:) = 0   !! NB: not used in the physics, only diagnostic...
667!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
668!!! for IDEALIZED CASES (NO NEED in pplay)
669pplev(subs,nlayer+1) = 0.
670!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
671
672! NOTE:
673! IF ( pplev(subs,nlayer+1) .le. 0 ) pplev(subs,nlayer+1)=ptop
674! cree des diagnostics delirants et aleatoires dans le transfert radiatif
675
676!---------!
677! Tracers !
678!---------!
679SELECT CASE (MARS_MODE)
680    CASE(0)  !! NO TRACERS (mars=0)
681    pq(subs,:,nq) = water_vapor_prof(:) !! NB: which is 0, actually (see above)
682    CASE(1)  !! WATER CYCLE (mars=1)
683    pq(subs,:,nq) = water_vapor_prof(:)
684    pq(subs,:,nq-1) = water_ice_prof(:)
685    CASE(2)  !! DUST CYCLE (mars=2)
686    pq(subs,:,nq) = water_vapor_prof(:) !! NB: which is 0, actually (see above)
687END SELECT
688
689ENDDO
690ENDDO
691
692!!*****************!!
693!! CLEAN THE PLACE !!
694!!*****************!!
695DEALLOCATE(q2_val)
696DEALLOCATE(qsurf_val)
697DEALLOCATE(tsoil_val)
698DEALLOCATE(dz8w_prof)
699DEALLOCATE(z_prof)
700DEALLOCATE(p8w_prof)
701DEALLOCATE(p_prof)
702DEALLOCATE(t_prof)
703DEALLOCATE(u_prof)
704DEALLOCATE(v_prof)
705DEALLOCATE(water_vapor_prof)
706DEALLOCATE(water_ice_prof)
707    !!! no use
708    !DEALLOCATE(pi_prof)
709    !DEALLOCATE(rho_prof)
710    !DEALLOCATE(th_prof)
711
712
713!!!!!!!!!!!!!!!!!!!!!!
714!!!!!!!!!!!!!!!!!!!!!!
715!! CALL LMD PHYSICS !! 
716!!!!!!!!!!!!!!!!!!!!!!
717!!!!!!!!!!!!!!!!!!!!!!
718
719
720!!***********************!!
721!! INIFIS, AT FIRST CALL !!
722!!***********************!!
723IF (firstcall .EQV. .true.) THEN
724print *, '** Mars ** LMD INITIALIZATION'
725include "../modif_mars/call_meso_inifis.inc"
726DEALLOCATE(aire_vec)
727DEALLOCATE(lat_vec)
728DEALLOCATE(lon_vec)
729DEALLOCATE(walbedodat)
730DEALLOCATE(winertiedat)
731DEALLOCATE(wphisfi)
732DEALLOCATE(wzmea)
733DEALLOCATE(wzstd)
734DEALLOCATE(wzsig)
735DEALLOCATE(wzgam)
736DEALLOCATE(wzthe)
737DEALLOCATE(wtheta)
738DEALLOCATE(wpsi)
739ENDIF
740        !! nearly obsolete
741        !print *, '** Mars ** Diagnostic files each ',wecri_phys,' phys. steps'
742        wecri_phys_sec=dt*float(wecri_phys)*float(wappel_phys)
743
744!!********!!
745!! PHYSIQ !!
746!!********!!
747call_physics : IF (wappel_phys .ne. 0) THEN
748
749!-------------------------------------------------------------------------------!
750! outputs:                                                                      !       
751!    pdu(ngrid,nlayermx)        \                                               !
752!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding    !
753!    pdt(ngrid,nlayermx)         /  variables due to physical processes.        !
754!    pdq(ngrid,nlayermx)        /                                               !
755!    pdpsrf(ngrid)             /                                                !
756!    tracerdyn                 call tracer in dynamical part of GCM ?           !
757!-------------------------------------------------------------------------------!
758print *, '** Mars ** CALL TO LMD PHYSICS'
759pdpsrf(:)=0.
760pdu(:,:)=0.
761pdv(:,:)=0.
762pdt(:,:)=0.
763pdq(:,:,:)=0.
764include "../modif_mars/call_meso_physiq.inc"
765DEALLOCATE(pplev)
766DEALLOCATE(pplay)
767DEALLOCATE(pphi)
768DEALLOCATE(pu)
769DEALLOCATE(pv)
770DEALLOCATE(pt)
771DEALLOCATE(pw)
772DEALLOCATE(pq)
773DEALLOCATE(wtsurf)
774DEALLOCATE(wco2ice)
775DEALLOCATE(wemis)
776DEALLOCATE(wq2)
777DEALLOCATE(wqsurf)
778DEALLOCATE(wtsoil)
779
780
781!-------------------------------!
782! PHYSIQ OUTPUT IN THE WRF FILE !
783!-------------------------------!
784DO j = jps,jpe
785DO i = ips,ipe
786subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
787#include "../modif_mars/module_lmd_driver_output3.inc" 
788       !  ^-- generated from Registry
789TSK(i,j) = output_tab2d(subs,ind_TSURF)
790ENDDO
791ENDDO
792DEALLOCATE(output_tab2d)
793DEALLOCATE(output_tab3d)
794
795!---------------------------------------------------------------------------------!
796! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
797!---------------------------------------------------------------------------------!
798dp_save(:)=pdpsrf(:)
799du_save(:,:)=pdu(:,:)
800dv_save(:,:)=pdv(:,:)
801dt_save(:,:)=pdt(:,:)
802dq_save(:,:,:)=pdq(:,:,:)
803
804ENDIF call_physics
805
806ENDIF    ! end of BIG LOOP 2.
807!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
808!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
809!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
810!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
811
812
813!!***************************!!
814!! DEDUCE TENDENCIES FOR WRF !!
815!!***************************!!
816RTHBLTEN(ims:ime,kms:kme,jms:jme)=0.
817RUBLTEN(ims:ime,kms:kme,jms:jme)=0.
818RVBLTEN(ims:ime,kms:kme,jms:jme)=0.
819PSFC(ims:ime,jms:jme)=p8w(ims:ime,kms,jms:jme) ! was done in surface driver in regular WRF
820!------------------------------------------------------------------!
821! WRF adds the physical and dynamical tendencies                   !
822! thus the tendencies must be passed as 'per second' tendencies    !
823! --when physics is not called, the applied physical tendency      !
824! --is the one calculated during the last call to physics          !
825!------------------------------------------------------------------!
826DO j = jps,jpe
827DO i = ips,ipe
828subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
829
830!------------!
831! zonal wind !
832!------------!
833RUBLTEN(i,kps:kpe,j) = pdu(subs,kps:kpe)                       
834
835!-----------------!
836! meridional wind !
837!-----------------!
838RVBLTEN(i,kps:kpe,j) = pdv(subs,kps:kpe)                       
839
840!-----------------------!
841! potential temperature !
842!-----------------------!
843! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
844RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe) / exner(i,kps:kpe,j)   
845
846!---------------------------!
847! update surface pressure   !
848! (cf CO2 cycle in physics) !
849!---------------------------!
850PSFC(i,j)=PSFC(i,j)+pdpsrf(subs) 
851
852!---------!
853! Tracers !
854!---------!
855SELECT CASE (MARS_MODE)
856CASE(0)
857        SCALAR(i,kps:kpe,j,:)=0.
858CASE(1)
859        !!! Water vapor
860        SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)
861        !!! Water ice
862        SCALAR(i,kps:kpe,j,3)=SCALAR(i,kps:kpe,j,3)+pdq(subs,kps:kpe,nq-1)
863CASE(2)
864        !!! Dust
865        SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)
866END SELECT
867
868        !!TODO: check if adding the whole tendency once, and set the
869        !!TODO: following tendencies to 0 until physics is called again
870        !!TODO: is a good strategy ?
871        !RUBLTEN(i,kps:kpe,j) = pdu(subs,kps:kpe)*ptimestep/dt
872        !RVBLTEN(i,kps:kpe,j) = pdv(subs,kps:kpe)*ptimestep/dt
873        !RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe)*ptimestep/dt
874        !RTHBLTEN(i,kps:kpe,j) = RTHBLTEN(i,kps:kpe,j)/exner(i,kps:kpe,j)
875        !PSFC(i,j)=PSFC(i,j)+pdpsrf(subs)*ptimestep/dt     
876        !SELECT CASE (MARS_MODE)
877        !CASE(0)
878        !SCALAR(i,kps:kpe,j,:)=0.
879        !CASE(1)
880        !!!! Water vapor
881        !SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*ptimestep/dt
882        !!!! Water ice
883        !SCALAR(i,kps:kpe,j,3)=SCALAR(i,kps:kpe,j,3)+pdq(subs,kps:kpe,nq-1)*ptimestep/dt
884        !CASE(2)
885        !!!! Dust
886        !SCALAR(i,kps:kpe,j,2)=SCALAR(i,kps:kpe,j,2)+pdq(subs,kps:kpe,nq)*ptimestep/dt
887        !END SELECT
888
889ENDDO
890ENDDO
891DEALLOCATE(pdpsrf)
892DEALLOCATE(pdu)
893DEALLOCATE(pdv)
894DEALLOCATE(pdt)
895DEALLOCATE(pdq)
896
897!!---------!
898!! display !
899!!---------!
900PRINT *, '** Mars ** Results from LMD physics'
901!PRINT *, 'u non-zero tendencies'
902!PRINT *, 'max',MAXVAL(RUBLTEN, MASK=RUBLTEN/=0.),&
903!        ' at',MAXLOC(RUBLTEN, MASK=RUBLTEN/=0.)
904!PRINT *, 'min',MINVAL(RUBLTEN, MASK=RUBLTEN/=0.),&
905!        ' at',MINLOC(RUBLTEN, MASK=RUBLTEN/=0.)
906!PRINT *, 'v non-zero tendencies'
907!PRINT *, 'max',MAXVAL(RVBLTEN, MASK=RVBLTEN/=0.),&
908!        ' at',MAXLOC(RVBLTEN, MASK=RVBLTEN/=0.)
909!PRINT *, 'min',MINVAL(RVBLTEN, MASK=RVBLTEN/=0.),&
910!        ' at',MINLOC(RVBLTEN, MASK=RVBLTEN/=0.)
911PRINT *, 'th non-zero tendencies'
912PRINT *, 'max',MAXVAL(RTHBLTEN, MASK=RTHBLTEN/=0.),&
913        ' at',MAXLOC(RTHBLTEN, MASK=RTHBLTEN/=0.)
914PRINT *, 'min',MINVAL(RTHBLTEN, MASK=RTHBLTEN/=0.),&
915        ' at',MINLOC(RTHBLTEN, MASK=RTHBLTEN/=0.)
916!!! STOP IF CRASH
917!IF (MAXVAL(RUBLTEN, MASK=RUBLTEN/=0.) == 0.) STOP
918!IF (MAXVAL(RVBLTEN, MASK=RVBLTEN/=0.) == 0.) STOP
919
920!!*****!!
921!! END !!
922!!*****!!
923print *,'** Mars ** END LMD PHYSICS'
924!----------------------------------------------------------------!
925! use debug (see solve_em) if you wanna display some message ... !
926!----------------------------------------------------------------!
927END SUBROUTINE lmd_driver
928
929!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
930      real function ls2sol(ls)
931
932!c  Returns solar longitude, Ls (in deg.), from day number (in sol),
933!c  where sol=0=Ls=0 at the northern hemisphere spring equinox
934
935      implicit none
936
937!c  Arguments:
938      real ls
939
940!c  Local:
941      double precision xref,zx0,zteta,zz
942!c      xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
943      double precision year_day
944      double precision peri_day,timeperi,e_elips
945      double precision pi,degrad
946      parameter (year_day=668.6d0) ! number of sols in a amartian year
947!c      data peri_day /485.0/
948      parameter (peri_day=485.35d0) ! date (in sols) of perihelion
949!c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
950      parameter (timeperi=1.90258341759902d0)
951      parameter (e_elips=0.0934d0)  ! eccentricity of orbit
952      parameter (pi=3.14159265358979d0)
953      parameter (degrad=57.2957795130823d0)
954
955      if (abs(ls).lt.1.0e-5) then
956         if (ls.ge.0.0) then
957            ls2sol = 0.0
958         else
959            ls2sol = year_day
960         end if
961         return
962      end if
963
964      zteta = ls/degrad + timeperi
965      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
966      xref = zx0-e_elips*dsin(zx0)
967      zz = xref/(2.*pi)
968      ls2sol = zz*year_day + peri_day
969      if (ls2sol.lt.0.0) ls2sol = ls2sol + year_day
970      if (ls2sol.ge.year_day) ls2sol = ls2sol - year_day
971
972      return
973      end function ls2sol
974!!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
975
976END MODULE module_lmd_driver
977
Note: See TracBrowser for help on using the repository browser.