source: trunk/WRF.COMMON/INTERFACES/dynphy_wrf_generic_lmd/update_inputs_physiq_mod.F

Last change on this file was 2074, checked in by aslmd, 6 years ago

correction in MESOSCALE generic update_inputs_physiq_mod, .eqv. should be used for booleans not .eq.
File size: 22.5 KB
Line 
1MODULE update_inputs_physiq_mod
2
3CONTAINS
4
5!SUBROUTINE update_inputs_physiq_time
6!SUBROUTINE update_inputs_physiq_tracers
7!SUBROUTINE update_inputs_physiq_constants
8!SUBROUTINE update_inputs_physiq_geom
9!SUBROUTINE update_inputs_physiq_surf
10!SUBROUTINE update_inputs_physiq_soil
11!SUBROUTINE update_inputs_physiq_turb
12!SUBROUTINE update_inputs_physiq_rad
13!SUBROUTINE update_inputs_physiq_slope
14
15!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
17SUBROUTINE update_inputs_physiq_time(&
18            JULYR,JULDAY,GMT,&
19            elaps,&
20            lct_input,lon_input,ls_input,&
21            MY)
22  USE variables_mod, only: JD_cur,JH_cur_split,phour_ini
23  use callkeys_mod, only : tlocked
24  INTEGER, INTENT(IN) :: JULDAY, JULYR
25  REAL, INTENT(IN) :: GMT,elaps,lon_input,ls_input,lct_input
26  REAL,INTENT(OUT) :: MY
27  REAL :: sec,nsec
28 
29  IF (JULYR .ne. 9999) THEN
30    if (tlocked .eqv. .false.) THEN
31      JH_cur_split = (GMT + elaps/3600.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT
32      JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type
33      JH_cur_split = JH_cur_split / 24.
34      JD_cur = (JULDAY - 1 + INT((3600.0*GMT+elaps)/86400))
35      JD_cur = MODULO(int(JD_cur),2)
36      MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600.0*GMT+elaps)/31968000
37      MY = INT(MY)
38    ELSE
39      JH_cur_split = (GMT)! + elaps/420000.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT
40      JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type
41      JH_cur_split = JH_cur_split / 24.
42      JD_cur = (JULDAY - 1 + INT((3600*GMT)/86400))!+elaps)/1.008e7))
43      JD_cur = MODULO(int(JD_cur),2)
44      MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600*GMT+elaps)/31968000
45      MY = INT(MY)
46    ENDIF
47  ELSE
48    if (tlocked .eqv. .false.) THEN
49      JH_cur_split = lct_input - lon_input / 15. + elaps/1500.0
50      JD_cur =  INT((sec*(lct_input - lon_input / 15.) + elaps)/36000)
51      !ptime = lct_input - lon_input / 15. + elaps/3600.
52      !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
53    ELSE
54      JH_cur_split = lct_input - lon_input / 15. !+ elaps/1500.0
55      !pday =  INT((sec*(lct_input - lon_input / 15.)+ elaps)/36000)
56      JD_cur =  INT((sec*(lct_input - lon_input / 15.))/3600)
57      !print*,'ptime',ptime
58      !print*,'pday',pday
59    !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
60      JH_cur_split = MODULO(ptime,24.)
61      JH_cur_split = JH_cur_split / 24.
62      JD_cur = MODULO(int(pday),365)
63      MY = 2024
64    ENDIF
65  ENDIF
66
67END SUBROUTINE update_inputs_physiq_time
68
69!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
70!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
71SUBROUTINE update_inputs_physiq_tracers(nq,MARS_MODE)
72
73  use tracer_h, only: noms,nqtot
74  !use phys_state_var_mod , only : qsurf
75  !USE dimphy, only : klon,klev
76  INTEGER, INTENT(IN) :: nq,MARS_MODE
77  INTEGER :: i,k
78!  CHARACTER(len=20), INTENT(INOUT) :: tname(nq) ! tracer names
79  logical :: end_of_file
80
81  !! TRACERS
82  IF (.not.ALLOCATED(noms)) ALLOCATE(noms(nq)) !! est fait dans initracer normalement
83     
84                 !! tableau dans tracer_mod.F90
85  print*,'nq',nq
86  print*,'MARS_MODE',MARS_MODE
87  nqtot=nq
88  if (MARS_MODE .eq. 1) THEN
89    noms(1)="h2o_vap"
90    noms(2)="h2o_ice"
91  else if (MARS_MODE .eq. 22)  then
92    noms(1)="co2_vap"
93    noms(2)="co2_ice"
94  else
95    noms(:)="zolbxs"
96  endif
97  !noms(:)=tname(:)
98  print*,'noms',noms
99  !!---------------------!!
100  !! OUTPUT FOR CHECKING !!
101  !!---------------------!!
102
103END SUBROUTINE update_inputs_physiq_tracers
104
105!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
106!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
107SUBROUTINE update_inputs_physiq_constants
108
109   !USE module_model_constants
110   !use comcstfi_h, only: omeg,mugaz
111   !use planete_h,  only: year_day,periheli,aphelie, &
112   !                       peri_day,obliquit,emin_turb, &
113   !                       lmixmin
114   use planete_mod, only: year_day, periastr, apoastr, peri_day,&
115                       obliquit, z0, lmixmin, emin_turb
116   use surfdat_h,  only: emissiv,iceradius, &
117                         emisice,dtemisice
118   !                       z0_default
119   !use comsoil_h, only: volcapa
120   use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r
121   !! comcstfi_h
122   use phys_state_var_mod, only :cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tsea_ice
123   !use iniorbit
124   !use time_phylmdz_mod, only: dtphys, daysec,day_ini
125 
126 
127     !open(17,file='controle.txt',form='formatted',status='old')
128     !rewind(17)
129     !read(17,*)
130     !read(17,*)
131     !read(17,*) day_ini !(tab0+3)
132     !read(17,*)
133     !read(17,*) !tab0+5)
134     !read(17,*) omeg !(tab0+6)
135     !read(17,*) !(tab0+7)
136     !read(17,*)  !(tab0+8)
137     !read(17,*)  !(tab0+9)
138     !read(17,*) daysec
139     !read(17,*) dtphys !tab0+11)
140     !read(17,*)
141     !read(17,*)
142     !read(17,*) year_day !(tab0+14)
143     !read(17,*) periastr !tab0+15)
144     !read(17,*) apoastr !tab0+16)
145     !read(17,*) peri_day !tab0+17)
146     !read(17,*) obliquit !tab0+18)
147     !read(17,*) z0
148     !read(17,*)
149     !read(17,*)
150     !read(17,*)
151     !read(17,*)
152     !read(17,*) emisice(1)
153     !read(17,*) emisice(2)
154     !read(17,*) emissiv
155     !read(17,*)
156     !read(17,*)
157     !read(17,*)
158     !read(17,*)
159     !read(17,*) iceradius(1)
160     !read(17,*) iceradius(2)
161     !read(17,*) dtemisice(1)
162     !read(17,*) dtemisice(2)
163     !close(17)
164     !cpp=(8.314511/(mugaz/1000.0))/rcp
165     !print*,'cpp',cpp
166     !print*,'g',g
167
168     !emissiv(:)=EMIS
169     !cloudfrac(:,:)=0.5
170     !totcloudfrac(:)=0.5
171     !hice(:)=0.
172     !rnat(:)=0.
173     !pctsrf_sic(:)=0.
174     !tsea_ice(:)=0.
175     !qsurf(:,:) = 0.
176     !print*,'iceradius',iceradius,'dtemisice',dtemisice
177     !print*,'apoastr,periastr,year_day,peri_day,obliq',apoastr,periastr,year_day,peri_day,obliquit
178     !print*,'emissiv',emissiv
179 
180     ! SUBROUTINE iniorbit(apoastr,periastr,year_day,peri_day,obliq)
181
182END SUBROUTINE update_inputs_physiq_constants
183
184!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
185!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186SUBROUTINE update_inputs_physiq_geom( &
187            ims,ime,jms,jme,&
188            ips,ipe,jps,jpe,&
189            JULYR,ngrid,nlayer,&
190            DX,DY,MSFT,&
191            lat_input, lon_input,&
192            XLAT,XLONG)
193
194   ! in WRF (share)
195   USE module_model_constants, only: DEGRAD,p0
196   ! in LMD (phymars)
197   !use comgeomfi_h, only: ini_fillgeom
198   ! in LMD (phy_common)
199   USE mod_grid_phy_lmdz, ONLY: init_grid_phy_lmdz
200   USE geometry_mod, ONLY: latitude,latitude_deg,&
201                           longitude,longitude_deg,&
202                           cell_area
203   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
204   !use planetwide_mod, only: planetwide_sumval
205   use comgeomfi_h, only: totarea, totarea_planet
206   use planete_mod, only: ini_planete_mod
207
208   INTEGER, INTENT(IN) :: ims,ime,jms,jme
209   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR,ngrid,nlayer
210   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
211     MSFT,XLAT,XLONG
212   REAL, INTENT(IN) :: dx,dy
213   REAL, INTENT(IN) :: lat_input, lon_input
214   INTEGER :: i,j,subs,ig,k
215   REAL,DIMENSION(ngrid) :: plon, plat, parea
216   REAL, DIMENSION(nlayer+1) :: znw
217   REAL*8, DIMENSION(nlayer+1) :: apdyn,bpdyn
218   REAL :: ptop
219   REAL*8 :: PP0
220   DO j = jps,jpe
221   DO i = ips,ipe
222
223    !-----------------------------------!
224    ! 1D subscript for physics "cursor" !
225    !-----------------------------------!
226    subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
227
228    !----------------------------------------!
229    ! Surface of each part of the grid (m^2) ! 
230    !----------------------------------------!
231    !parea(subs) = dx*dy                           !! 1. idealized cases - computational grid
232    parea(subs) = (dx/msft(i,j))*(dy/msft(i,j))    !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
233    !parea(subs)=dx*dy/msfu(i,j)                   !! 3. special for Mercator GCM-like simulations
234
235    !---------------------------------------------!
236    ! Mass-point latitude and longitude (radians) !
237    !---------------------------------------------!
238    IF (JULYR .ne. 9999) THEN
239     plat(subs) = XLAT(i,j)*DEGRAD
240     plon(subs) = XLONG(i,j)*DEGRAD
241    ELSE
242     !!! IDEALIZED CASE
243     IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lat: ',lat_input
244     IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lon: ',lon_input
245     plat(subs) = lat_input*DEGRAD
246     plon(subs) = lon_input*DEGRAD
247    ENDIF
248
249   ENDDO
250   ENDDO
251
252   !! FILL GEOMETRICAL ARRAYS !!
253   !call ini_fillgeom(ngrid,plat,plon,parea)
254
255   !!! ----------------------------------------------------------
256   !!! --- initializing geometry in phy_common
257   !!! --- (this is quite planet-independent)
258   !!! ----------------------------------------------------------
259   ! initialize mod_grid_phy_lmdz
260   CALL init_grid_phy_lmdz(1,1,ipe-ips+1,jpe-jps+1,nlayer)
261   ! fill in geometry_mod variables
262   ! ... copy over local grid longitudes and latitudes
263   ! ... partly what is done in init_geometry
264   IF(.not.ALLOCATED(longitude)) ALLOCATE(longitude(ngrid))
265   IF(.not.ALLOCATED(longitude_deg)) ALLOCATE(longitude_deg(ngrid))
266   IF(.not.ALLOCATED(latitude)) ALLOCATE(latitude(ngrid))
267   IF(.not.ALLOCATED(latitude_deg)) ALLOCATE(latitude_deg(ngrid))
268   IF(.not.ALLOCATED(cell_area)) ALLOCATE(cell_area(ngrid))
269   longitude(:) = plon(:)
270   latitude(:) = plat(:)
271   longitude_deg(:) = plon(:)/DEGRAD
272   latitude_deg(:) = plat(:)/DEGRAD
273   cell_area(:) = parea(:)
274   !call planetwide_sumval(parea,totarea_planet)
275   !print*,'parea',parea(1)
276   !totarea=SSUM(ngrid,parea,1)
277   totarea=ngrid*parea(1)
278   !totarea_planet=SSUM(ngrid,parea,1)
279   totarea_planet=ngrid*parea(1)
280
281   IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
282   IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
283   IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
284   IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
285   DO ig=1,ngrid
286     sinlat(ig)=sin(plat(ig))
287     coslat(ig)=cos(plat(ig))
288     sinlon(ig)=sin(plon(ig))
289     coslon(ig)=cos(plon(ig))
290   ENDDO
291   
292   open(unit=12,file='levels',form='formatted',status='old')
293   rewind(12)
294   DO k=1, nlayer
295   read(12,*) znw(k)
296   !write(6,*) 'read level ', k,grid%znw(k)
297   ENDDO
298   close(12)
299   
300   ptop=0.5
301   apdyn=ptop*(1-znw)
302   bpdyn=znw
303   PP0=p0
304   CALL ini_planete_mod(nlayer,PP0,apdyn,bpdyn)
305   !!! ----------------------------------------------------------
306
307END SUBROUTINE update_inputs_physiq_geom
308
309!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
311SUBROUTINE update_inputs_physiq_surf( &
312            ims,ime,jms,jme,&
313            ips,ipe,jps,jpe,&
314            JULYR,MARS_MODE,&
315            M_ALBEDO,CST_AL,&
316            M_TSURF,M_EMISS,M_CO2ICE,&
317            M_GW,M_Z0,CST_Z0,&
318            M_H2OICE,&
319            phisfi_val)
320
321   use surfdat_h, only: phisfi, albedodat,  &
322                       zmea, zstd, zsig, zgam, zthe
323   use slab_ice_h, only: noceanmx
324   use planete_mod, only: z0
325   use phys_state_var_mod, only : tsurf, emis, qsurf,tslab
326
327   INTEGER, INTENT(IN) :: ims,ime,jms,jme
328   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR,MARS_MODE
329   INTEGER :: i,j,subs,nlast,iq
330   REAL, INTENT(IN  ) :: CST_AL, phisfi_val, CST_Z0
331   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
332     M_ALBEDO,M_TSURF,M_EMISS,M_CO2ICE,M_H2OICE,M_Z0
333   REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: M_GW 
334
335   !print*,'ALLOCATED(phisfi)',ALLOCATED(phisfi)
336   !print*,'size phisfi',size(phisfi)
337   DO j = jps,jpe
338   DO i = ips,ipe
339
340     !-----------------------------------!
341     ! 1D subscript for physics "cursor" !
342     !-----------------------------------!
343     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
344
345     !---------------------!
346     ! Ground geopotential !
347     !---------------------!
348     phisfi(subs) = phisfi_val
349!     print*,'size phisfi',size(phisfi)
350     !print*,'phisfi',phisfi(subs)
351     !---------------!
352     ! Ground albedo !
353     !---------------!
354     IF (JULYR .ne. 9999) THEN
355      IF (CST_AL == 0) THEN
356       albedodat(subs)=M_ALBEDO(i,j)
357      ELSE
358       albedodat(subs)=CST_AL
359       IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT ALBEDO ', albedodat
360      ENDIF
361     ELSE
362      IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION albedo: ', CST_AL
363      albedodat(subs)=CST_AL
364     ENDIF
365
366     !-----------------------------------------!
367     ! Gravity wave parametrization            !
368     ! NB: usually 0 in mesoscale applications !
369     !-----------------------------------------!
370     zmea(subs)=0.
371     zstd(subs)=0.
372     zsig(subs)=0.
373     zgam(subs)=0.
374     zthe(subs)=0.
375
376     !----------------------------!
377     ! Variable surface roughness !
378     !----------------------------!
379     z0 = CST_Z0
380     !IF (JULYR .ne. 9999) THEN
381     ! IF (CST_Z0 == 0) THEN
382     !  z0(subs) = M_Z0(i,j)
383     ! ELSE
384     !  z0(subs) = CST_Z0
385     !  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT SURF ROUGHNESS (m) ',CST_Z0
386     ! ENDIF
387     !ELSE
388     ! IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION z0 (m) ', CST_Z0
389     ! z0(subs)=CST_Z0
390     !ENDIF
391     !!!!! ADDITIONAL SECURITY. THIS MIGHT HAPPEN WITH OLD INIT FILES.
392     !IF (z0(subs) == 0.) THEN
393     !  IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
394     !  z0(subs) = 0.01
395     !ENDIF
396     !!!!! ADDITIONAL SECURITY. INTERP+SMOOTH IN GEOGRID MIGHT YIELD NEGATIVE Z0 !!!
397     !IF (z0(subs) < 0.) THEN
398     !  PRINT *, 'WELL, z0 is NEGATIVE, this is impossible. better stop here.'
399     !  PRINT *, 'advice --> correct interpolation / smoothing of z0 in WPS'
400     !  PRINT *, '           -- or check the constant value set in namelist.input'
401     !  STOP
402     !ENDIF
403
404     !-----------------------------------------------!
405     ! Ground temperature, emissivity, CO2 ice cover !
406     !-----------------------------------------------!
407     tsurf(subs) = M_TSURF(i,j)
408     emis(subs) = M_EMISS(i,j)     
409     !do i=1,noceanmx
410     tslab(subs,:)=tsurf(subs)
411     !enddo
412     !-------------------!
413     ! Tracer at surface !
414     !-------------------!
415     qsurf(subs,:)=0. ! default case
416     SELECT CASE (MARS_MODE)
417       CASE(1)
418       qsurf(subs,2)=M_H2OICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
419                                      !! ----- retrocompatible ancienne physique
420                                      !! ----- [H2O ice is last tracer in qsurf in LMD physics]
421       CASE(2)   
422       qsurf(subs,1)=0.                !! not coupled with lifting for the moment [non remobilise]
423       !CASE(3)
424       !qsurf(subs,1)=q_prof(1,1)                !!! temporaire, a definir       
425       !qsurf(subs,2)=q_prof(1,2)     
426       CASE(11)
427       qsurf(subs,2)=M_H2OICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
428       qsurf(subs,3)=0.                !! not coupled with lifting for the moment [non remobilise]
429       CASE(12)
430       qsurf(subs,2)=M_H2OICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
431       qsurf(subs,3)=0.                !! not coupled with lifting for the moment [non remobilise]
432     END SELECT
433
434   ENDDO
435   ENDDO
436 
437   !!---------------------!!
438   !! OUTPUT FOR CHECKING !!
439   !!---------------------!!
440   nlast = (ipe-ips+1)*(jpe-jps+1)
441   print*,"check: phisfi",phisfi(1),phisfi(nlast)
442   print*,"check: albedodat",albedodat(1),albedodat(nlast)
443   print*,"check: zmea",zmea(1),zmea(nlast)
444   print*,"check: zstd",zstd(1),zstd(nlast)
445   print*,"check: zsig",zsig(1),zsig(nlast)
446   print*,"check: zgam",zgam(1),zgam(nlast)
447   print*,"check: zthe",zthe(1),zthe(nlast)
448   print*,"check: z0",z0
449   print*,"check: tsurf",tsurf(1),tsurf(nlast)
450   print*,"check: emis",emis(1),emis(nlast)
451   print*,"check: qsurf",qsurf(1,:),qsurf(nlast,:)
452
453END SUBROUTINE update_inputs_physiq_surf
454
455!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
457SUBROUTINE update_inputs_physiq_soil( &
458            ims,ime,jms,jme,&
459            ips,ipe,jps,jpe,&
460            JULYR,nsoil,&
461            M_TI,CST_TI,&
462            M_ISOIL,M_DSOIL,&
463            M_TSOIL,M_TSURF)
464
465   use comsoil_h, only: inertiedat,mlayer,layer,volcapa
466   use phys_state_var_mod, only: tsoil
467
468   INTEGER, INTENT(IN) :: ims,ime,jms,jme
469   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR,nsoil
470   INTEGER :: i,j,subs,nlast
471   REAL, INTENT(IN  ) :: CST_TI
472   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
473     M_TI, M_TSURF
474   REAL, DIMENSION( ims:ime, nsoil, jms:jme ), INTENT(IN)  :: &
475     M_TSOIL, M_ISOIL, M_DSOIL
476   REAL :: inertiedat_val
477   REAL :: lay1,alpha
478
479   DO j = jps,jpe
480   DO i = ips,ipe
481
482     !-----------------------------------!
483     ! 1D subscript for physics "cursor" !
484     !-----------------------------------!
485     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
486
487     !-----------------!
488     ! Thermal Inertia !
489     !-----------------!
490     IF (JULYR .ne. 9999) THEN
491      IF (CST_TI == 0) THEN
492       inertiedat_val=M_TI(i,j)
493      ELSE
494       inertiedat_val=CST_TI
495       IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT THERMAL INERTIA ', inertiedat_val
496      ENDIF
497     ELSE
498      IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION inertia: ',CST_TI
499      inertiedat_val=CST_TI
500     ENDIF
501     !inertiedat(subs) = inertiedat_val
502     !--pb de dimensions???!!???
503     IF (JULYR .ne. 9999) THEN
504       inertiedat(subs,:)=M_ISOIL(i,:,j) !! verifier que cest bien hires TI en surface
505       mlayer(0:nsoil-1)=M_DSOIL(i,:,j)
506     ELSE
507        IF ( nsoil .lt. 18 ) THEN
508            PRINT *,'** Mars ** WRONG NUMBER OF SOIL LAYERS. SHOULD BE 18 AND IT IS ',nsoil
509            STOP
510        ENDIF
511        IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION isoil and dsoil standard'
512        do k=1,nsoil
513         inertiedat(subs,k) = inertiedat_val
514         !mlayer(k-1) = sqrt(887.75/3.14)*((2.**(k-0.5))-1.) * inertiedat_val / wvolcapa    ! old setting
515         mlayer(k-1) = 2.E-4 * (2.**(k-0.5-1.))                                            ! new gcm settings
516        enddo
517     ENDIF
518     IF ( (i == ips) .AND. (j == jps) ) &
519         PRINT *,'** Mars ** TI and depth profiles are',inertiedat(subs,:)!,mlayer(0:nsoil-1)
520
521     !!!!!!!!!!!!!!!!! DONE in soil_setting.F
522     ! 1.5 Build layer(); following the same law as mlayer()
523     ! Assuming layer distribution follows mid-layer law:
524     ! layer(k)=lay1*alpha**(k-1)
525     lay1=sqrt(mlayer(0)*mlayer(1))
526     alpha=mlayer(1)/mlayer(0)
527     do k=1,nsoil
528       layer(k)=lay1*(alpha**(k-1))
529     enddo
530
531     !------------------------!
532     ! Deep soil temperatures !
533     !------------------------!
534     IF (M_TSOIL(i,1,j) .gt. 0. .and. JULYR .ne. 9999) THEN
535       tsoil(subs,:)=M_TSOIL(i,:,j)
536     ELSE
537       IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** no tsoil. set it to tsurf.'
538       do k=1,nsoil
539        !print*,'M_TSURF(i,j)',M_TSURF(1,:)
540        !print*,'size M_TSURF',size(M_TSURF)
541        !print*,'size tsoil',size(tsoil)
542        tsoil(subs,k) = M_TSURF(i,j)
543        !print*,'tsoil(subs,k)',tsoil(subs,k)
544       enddo
545     ENDIF
546
547   ENDDO
548   ENDDO
549
550   volcapa=1.e6   
551
552   print*,'zolbxs'
553   !!---------------------!!
554   !! OUTPUT FOR CHECKING !!
555   !!---------------------!!
556   nlast = (ipe-ips+1)*(jpe-jps+1)
557   print*,"check: inertiedat",inertiedat(1,:),inertiedat(nlast,:)
558   print*,"check: mlayer",mlayer(:)
559   print*,"check: layer",layer(:)
560   print*,"check: tsoil",tsoil(1,:),tsoil(nlast,:)
561
562END SUBROUTINE update_inputs_physiq_soil
563
564!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
565!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
566SUBROUTINE update_inputs_physiq_turb( &
567            ims,ime,jms,jme,kms,kme,&
568            ips,ipe,jps,jpe,&
569            RESTART,isles,&
570            M_Q2,M_WSTAR)
571
572   use turb_mod, only: q2,wstar,turb_resolved
573   !use phys_state_var_mod, only : q2,
574
575   INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme
576   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe
577   INTEGER :: i,j,subs,nlast     
578   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: M_WSTAR
579   REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(IN) :: M_Q2
580   LOGICAL, INTENT(IN ) :: RESTART,isles
581
582   !! to know if this is turbulence-resolving run or not
583
584   turb_resolved = isles
585   print*,'isles',isles
586
587   DO j = jps,jpe
588   DO i = ips,ipe
589
590     !-----------------------------------!
591     ! 1D subscript for physics "cursor" !
592     !-----------------------------------!
593     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
594
595     !PBL wind variance
596     IF (.not. restart) THEN
597      q2(subs,:) = 1.E-6     
598      wstar(subs)=0.
599     ELSE
600      q2(subs,:)=M_Q2(i,:,j)!
601      !q2(subs,:) = 1.e-3
602      wstar(subs)=M_WSTAR(i,j)
603     ENDIF
604
605   ENDDO
606   ENDDO
607 
608   !!---------------------!!
609   !! OUTPUT FOR CHECKING !!
610   !!---------------------!!
611   nlast = (ipe-ips+1)*(jpe-jps+1)
612   print*,"check: q2",q2(1,1)!,q2(nlast,:)
613   print*,"check: wstar",wstar(1),wstar(nlast)
614
615END SUBROUTINE update_inputs_physiq_turb
616
617!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
618!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
619SUBROUTINE update_inputs_physiq_rad( &
620            ims,ime,jms,jme,&
621            ips,ipe,jps,jpe,&
622            RESTART,&
623            M_FLUXRAD)
624
625   !use dimradmars_mod, only: fluxrad
626   use phys_state_var_mod, only : fluxrad
627
628   INTEGER, INTENT(IN) :: ims,ime,jms,jme
629   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe
630   INTEGER :: i,j,subs,nlast     
631   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: M_FLUXRAD
632   LOGICAL, INTENT(IN ) :: RESTART
633
634   DO j = jps,jpe
635   DO i = ips,ipe
636
637     !-----------------------------------!
638     ! 1D subscript for physics "cursor" !
639     !-----------------------------------!
640     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
641
642     ! fluxrad_save
643     IF (.not. restart) THEN
644      fluxrad(subs)=0.
645     ELSE
646      fluxrad(subs)=M_FLUXRAD(i,j)
647     ENDIF
648     !! et fluxrad_sky ???!???
649
650   ENDDO
651   ENDDO
652 
653   !!---------------------!!
654   !! OUTPUT FOR CHECKING !!
655   !!---------------------!!
656   nlast = (ipe-ips+1)*(jpe-jps+1)
657   print*,"check: fluxrad",fluxrad(1),fluxrad(nlast)
658
659END SUBROUTINE update_inputs_physiq_rad
660
661!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
662!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
663SUBROUTINE update_inputs_physiq_slope( &
664            ims,ime,jms,jme,&
665            ips,ipe,jps,jpe,&
666            JULYR,&
667            SLPX,SLPY)
668
669   !USE module_model_constants, only: DEGRAD
670   !USE slope_mod, ONLY: theta_sl, psi_sl
671
672   INTEGER, INTENT(IN) :: ims,ime,jms,jme
673   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR
674   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: SLPX,SLPY
675   INTEGER :: i,j,subs,nlast
676
677END SUBROUTINE update_inputs_physiq_slope
678
679END MODULE update_inputs_physiq_mod
Note: See TracBrowser for help on using the repository browser.