source: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 @ 2903

Last change on this file since 2903 was 2896, checked in by romain.vande, 2 years ago

Mars PCM:
First modifications to introduce subslope parametrization in surface processes.
New variables are added in the new module comslope_mod to describe the subslopes:

_ nslope (number of subslope)
_ subslope_dist (distribution of the slopes)
_ def_slope (bound of the slope statistics / repartitions) etc..

These variables are added to the starfi file.
Other GCM variables, outputs and subroutines are not modified yet.
Only nslope=1 is accepted so far (ie. same as no subslope parametrization)
File size: 25.0 KB
Line 
1module phyetat0_mod
2
3implicit none
4  real,save :: tab_cntrl_mod(100)
5
6!$OMP THREADPRIVATE(tab_cntrl_mod)
7
8contains
9
10subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
11                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf, &
12                     tauscaling,totcloudfrac,wstar,watercap,def_slope, &
13                     def_slope_mean,subslope_dist)
14
15  use tracer_mod, only: noms ! tracer names
16  use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
17                       zmea, zstd, zsig, zgam, zthe, hmons, summit, base,&
18                       watercaptag
19  use iostart, only: nid_start, open_startphy, close_startphy, &
20                     get_field, get_var, inquire_field, &
21                     inquire_dimension, inquire_dimension_length
22  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
23  use compute_dtau_mod, only: dtau
24  use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
25  use dust_param_mod, only: dustscaling_mode
26  USE ioipsl_getin_p_mod, ONLY : getin_p
27  use comsoil_h, only: flux_geo
28  USE comslope_mod, ONLY: nslope, major_slope
29  implicit none
30 
31  include "callkeys.h"
32!======================================================================
33! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
34!  Adaptation � Mars : Yann Wanherdrick
35! Objet: Lecture de l etat initial pour la physique
36! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
37!                      r4 or r8 restarts independently of having compiled
38!                      the GCM in r4 or r8)
39!         June 2013 TN : Possibility to read files with a time axis
40!         November 2013 EM : Enabeling parallel, using iostart module
41!         March 2020 AD: Enabling initialization of physics without startfi
42!                        flag: startphy_file
43!======================================================================
44  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
45  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
46!======================================================================
47!  Arguments:
48!  ---------
49!  inputs:
50!  logical,intent(in) :: startphy_file ! .true. if reading start file
51  character*(*),intent(in) :: fichnom ! "startfi.nc" file
52  integer,intent(in) :: tab0
53  integer,intent(in) :: Lmodif
54  integer,intent(in) :: nsoil ! # of soil layers
55  integer,intent(in) :: ngrid ! # of atmospheric columns
56  integer,intent(in) :: nlay ! # of atmospheric layers
57  integer,intent(in) :: nq
58  integer :: day_ini
59  real :: time0
60
61!  outputs:
62  real,intent(out) :: tsurf(ngrid) ! surface temperature
63  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
64  real,intent(out) :: albedo(ngrid,2) ! surface albedo
65  real,intent(out) :: emis(ngrid) ! surface emissivity
66  real,intent(out) :: q2(ngrid,nlay+1) !
67  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
68  real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
69  real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
70  real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
71  real,intent(out) :: watercap(ngrid) ! h2o_ice_cover
72  real, intent(out) :: def_slope(nslope+1) !boundaries for bining of the slopes
73  real, intent(out) :: def_slope_mean(nslope)
74  real, intent(out) :: subslope_dist(ngrid,nslope) !undermesh statistics
75!======================================================================
76!  Local variables:
77
78      real surffield(ngrid) ! to temporarily store a surface field
79      real xmin,xmax ! to display min and max of a field
80!
81      INTEGER ig,iq,lmax,islope
82      INTEGER nid, nvarid
83      INTEGER ierr, i, nsrf
84!      integer isoil
85!      INTEGER length
86!      PARAMETER (length=100)
87      CHARACTER*7 str7
88      CHARACTER*2 str2
89      CHARACTER*1 yes
90!
91      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
92      INTEGER nqold
93
94! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
95      logical :: oldtracernames=.false.
96      integer :: count
97      character(len=30) :: txt ! to store some text
98
99! specific for time
100      REAL,ALLOCATABLE :: time(:) ! times stored in start
101      INTEGER timelen ! number of times stored in the file
102      INTEGER indextime ! index of selected time
103
104      INTEGER :: edges(3),corner(3)
105      LOGICAL :: found
106
107      REAL :: timestart ! to pick which initial state to start from
108      REAL :: surfemis  ! constant emissivity when no startfi
109      REAL :: surfalbedo  ! constant albedo when no startfi
110
111      REAL :: watercaptag_tmp(ngrid)
112
113! Sub-grid scale slopes
114      LOGICAL :: startphy_slope !to be retrocompatible and add the nslope dimension
115      REAL, ALLOCATABLE :: default_def_slope(:)
116      REAL :: sum_dist
117      REAL :: current_max !var to find max distrib slope
118
119      CHARACTER(len=5) :: modname="phyetat0"
120
121write(*,*) "phyetat0: startphy_file", startphy_file
122
123if (startphy_file) then
124   ! open physics initial state file:
125   call open_startphy(fichnom)
126   ! possibility to modify tab_cntrl in tabfi
127   write(*,*)
128   write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
129   call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
130               p_omeg,p_g,p_mugaz,p_daysec,time0)
131else ! "academic" initialization of planetary parameters via tabfi
132   call tabfi (0,0,0,day_ini,lmax,p_rad, &
133               p_omeg,p_g,p_mugaz,p_daysec,time0)
134endif ! of if (startphy_file)
135
136if(nslope.ne.1) then
137  call abort_physic(modname, &
138                "phyetat0: For now, nslope should be 1 (set in comslope_mod)",1)
139endif
140
141allocate(default_def_slope(nslope+1))
142!Sub-grid scale  subslopes
143if (nslope.eq.7) then
144  default_def_slope(1) = -43.
145  default_def_slope(2) = -19.
146  default_def_slope(3) = -9.
147  default_def_slope(4) = -3.
148  default_def_slope(5) = 3.
149  default_def_slope(6) = 9.
150  default_def_slope(7) = 19.
151  default_def_slope(8) = 43.
152elseif (nslope.eq.5) then
153  default_def_slope(1) = -43.
154  default_def_slope(2) = -9.
155  default_def_slope(3) = -3.
156  default_def_slope(4) = 3.
157  default_def_slope(5) = 9.
158  default_def_slope(6) = 43.
159elseif (nslope.eq.1) then
160  default_def_slope(1) = 0.
161  default_def_slope(2) = 0.
162endif
163
164if (startphy_file) then
165 call get_var("def_slope",def_slope,found)
166   if(.not.found) then
167     startphy_slope=.false.
168     write(*,*)'slope_settings: Problem while reading <def_slope>'
169     write(*,*)'default def_slope will be used'
170     do islope=1,nslope+1
171       def_slope(islope) = default_def_slope(islope)
172     enddo
173     write(*,*)'computing corresponding distribution <subslope_dist>'
174     write(*,*)'For now, woth nslope=1, subslope_dist is straigforward'
175     write(*,*)'Later this operation will be done by newstart with a specific routine'
176     subslope_dist(:,:)=1.
177     !call subslope_mola(ngrid,nslope,def_slope,subslope_dist)
178   else
179     startphy_slope=.true.
180     call get_field("subslope_dist",subslope_dist,found,indextime)
181     if(.not.found) then
182       write(*,*)'slope_settings: Problem while reading <subslope_dist>'
183       write(*,*)'computing a new distribution'
184       write(*,*)'For now, woth nslope=1, subslope_dist is straigforward'
185       write(*,*)'Later this operation will be done by newstart with a specific routine'
186       subslope_dist(:,:)=1.
187     !call subslope_mola(ngrid,nslope,def_slope,subslope_dist)
188     endif
189    endif
190else ! startphy_file 
191 do islope=1,nslope+1
192       def_slope(islope) = default_def_slope(islope)
193     enddo
194     write(*,*)'computing corresponding distribution <subslope_dist>'
195     write(*,*)'For now, woth nslope=1, subslope_dist is straigforward'
196     write(*,*)'Later this operation will be done by newstart with a specific routine'
197     subslope_dist(:,:)=1.
198     !call subslope_mola(ngrid,nslope,def_slope,subslope_dist)
199endif   
200
201do islope=1,nslope
202    def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
203enddo
204 
205DO ig = 1,ngrid
206  sum_dist = 0.
207  DO islope = 1,nslope
208     sum_dist = sum_dist + subslope_dist(ig,islope)
209  ENDDO
210  DO islope = 1,nslope
211    subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
212  ENDDO
213ENDDO
214
215!Now determine the major subslope, ie. the maximal distribution
216
217DO ig=1,ngrid
218  major_slope(ig)=1
219  current_max=subslope_dist(ig,1)
220  DO islope=2,nslope
221    if(subslope_dist(ig,islope).GT.current_max) then
222      major_slope(ig)=islope
223      current_max=subslope_dist(ig,islope)
224    ENDIF
225  ENDDO
226ENDDO
227
228if (startphy_file) then
229   ! Load surface geopotential:
230   call get_field("phisfi",phisfi,found)
231   if (.not.found) then
232     call abort_physic(modname, &
233                "phyetat0: Failed loading <phisfi>",1)
234   endif
235else
236  phisfi(:)=0.
237endif ! of if (startphy_file)
238write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
239               minval(phisfi), maxval(phisfi)
240
241
242if (startphy_file) then
243   ! Load bare ground albedo:
244   call get_field("albedodat",albedodat,found)
245   if (.not.found) then
246     call abort_physic(modname, &
247                "phyetat0: Failed loading <albedodat>",1)
248   endif
249else ! If no startfi file, use parameter surfalbedo in def file
250  surfalbedo=0.1
251  call getin_p("surfalbedo_without_startfi",surfalbedo)
252  print*,"surfalbedo_without_startfi",surfalbedo
253  albedodat(:)=surfalbedo
254endif ! of if (startphy_file)
255write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
256             minval(albedodat), maxval(albedodat)
257
258! ZMEA
259if (startphy_file) then
260   call get_field("ZMEA",zmea,found)
261   if (.not.found) then
262     call abort_physic(modname, &
263                "phyetat0: Failed loading <ZMEA>",1)
264   endif
265else
266  zmea(:)=0.
267endif ! of if (startphy_file)
268write(*,*) "phyetat0: <ZMEA> range:", &
269            minval(zmea), maxval(zmea)
270
271
272! ZSTD
273if (startphy_file) then
274   call get_field("ZSTD",zstd,found)
275   if (.not.found) then
276     call abort_physic(modname, &
277                "phyetat0: Failed loading <ZSTD>",1)
278   endif
279else
280  zstd(:)=0.
281endif ! of if (startphy_file)
282write(*,*) "phyetat0: <ZSTD> range:", &
283            minval(zstd), maxval(zstd)
284
285
286! ZSIG
287if (startphy_file) then
288   call get_field("ZSIG",zsig,found)
289   if (.not.found) then
290     call abort_physic(modname, &
291                "phyetat0: Failed loading <ZSIG>",1)
292   endif
293else
294  zsig(:)=0.
295endif ! of if (startphy_file)
296write(*,*) "phyetat0: <ZSIG> range:", &
297            minval(zsig), maxval(zsig)
298
299
300! ZGAM
301if (startphy_file) then
302   call get_field("ZGAM",zgam,found)
303   if (.not.found) then
304     call abort_physic(modname, &
305                "phyetat0: Failed loading <ZGAM>",1)
306   endif
307else
308  zgam(:)=0.
309endif ! of if (startphy_file)
310write(*,*) "phyetat0: <ZGAM> range:", &
311            minval(zgam), maxval(zgam)
312
313
314! ZTHE
315if (startphy_file) then
316   call get_field("ZTHE",zthe,found)
317   if (.not.found) then
318     call abort_physic(modname, &
319                "phyetat0: Failed loading <ZTHE>",1)
320   endif
321else
322  zthe(:)=0.
323endif ! of if (startphy_file)
324write(*,*) "phyetat0: <ZTHE> range:", &
325             minval(zthe), maxval(zthe)
326
327! hmons
328if (startphy_file) then
329   call get_field("hmons",hmons,found)
330   if (.not.found) then
331     write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
332     if (rdstorm) then
333     call abort_physic(modname, &
334                "phyetat0: Failed loading <hmons>",1)
335     else
336       write(*,*) "will continue anyway..."
337       write(*,*) "because you may not need it."
338       hmons(:)=0.
339     end if ! if (rdstorm)
340   else
341     do ig=1,ngrid
342       if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
343     enddo
344   endif ! (.not.found)
345else
346   hmons(:)=0.
347endif ! if (startphy_file)
348write(*,*) "phyetat0: <hmons> range:", &
349            minval(hmons), maxval(hmons)
350
351
352! summit
353if (startphy_file) then
354   call get_field("summit",summit,found)
355   if (.not.found) then
356     write(*,*) "WARNING: phyetat0: Failed loading <summit>"
357     if (rdstorm) then
358     call abort_physic(modname, &
359                "phyetat0: Failed loading <summit>",1)
360     else
361       write(*,*) "will continue anyway..."
362       write(*,*) "because you may not need it."
363       summit(:)=0.
364     end if
365   else
366     do ig=1,ngrid
367       if (summit(ig) .eq. -999999.)  summit(ig)=0.
368     enddo
369   endif ! if (.not.found)
370else
371   summit(:)=0. 
372endif ! if (startphy_file)
373write(*,*) "phyetat0: <summit> range:", &
374            minval(summit), maxval(summit)
375
376
377! base
378if (startphy_file) then
379   call get_field("base",base,found)
380   if (.not.found) then
381     write(*,*) "WARNING: phyetat0: Failed loading <base>"
382     if (rdstorm) then
383     call abort_physic(modname, &
384                "phyetat0: Failed loading <base>",1)
385     else
386       write(*,*) "will continue anyway..."
387       write(*,*) "because you may not need it."
388       base(:)=0.
389     end if
390   else
391     do ig=1,ngrid
392       if (base(ig) .eq. -999999.)  base(ig)=0.
393     enddo
394   endif ! if(.not.found)
395else
396   base(:)=0.
397endif ! if (startphy_file)
398write(*,*) "phyetat0: <base> range:", &
399            minval(base), maxval(base)
400
401! Time axis
402! obtain timestart from run.def
403timestart=-9999 ! default value
404call getin_p("timestart",timestart)
405if (startphy_file) then
406   found=inquire_dimension("Time")
407   if (.not.found) then
408     indextime = 1
409     write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
410   else
411     write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
412     timelen=inquire_dimension_length("Time")
413     allocate(time(timelen))
414     ! load "Time" array:
415     call get_var("Time",time,found)
416     if (.not.found) then
417     call abort_physic(modname, &
418                "phyetat0: Failed loading <Time>",1)
419     endif
420     ! seclect the desired time index
421     IF (timestart .lt. 0) THEN  ! default: we use the last time value
422       indextime = timelen
423     ELSE  ! else we look for the desired value in the time axis
424       indextime = 0
425       DO i=1,timelen
426         IF (abs(time(i) - timestart) .lt. 0.01) THEN
427           indextime = i
428           EXIT
429         ENDIF
430       ENDDO
431       IF (indextime .eq. 0) THEN
432         PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
433         PRINT*, "Stored times are:"
434         DO i=1,timelen
435            PRINT*, time(i)
436         ENDDO
437         call abort_physic(modname,"phyetat0: Time error",1)
438       ENDIF
439     ENDIF ! of IF (timestart .lt. 0)
440     ! In startfi the absolute date is day_ini + time0 + time
441     ! For now on, in the GCM physics, it is day_ini + time0
442     time0 = time(indextime) + time0
443     day_ini = day_ini + INT(time0)
444     time0 = time0 - INT(time0)   
445     PRINT*, "phyetat0: Selected time ",time(indextime), &
446             " at index ",indextime
447     DEALLOCATE(time)
448   endif ! of if Time not found in file
449
450   call ini_tab_controle_dyn_xios(day_ini)
451
452else
453  indextime = 1
454endif ! if (startphy_file)
455
456! Dust conversion factor
457if (startphy_file) then
458   call get_field("tauscaling",tauscaling,found,indextime)
459   if (.not.found) then
460     write(*,*) "phyetat0: <tauscaling> not in file"
461     tauscaling(:) = 1
462   endif
463else
464   tauscaling(:) = 1
465endif ! if (startphy_file)
466write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
467            minval(tauscaling), maxval(tauscaling)
468
469! dust_rad_adjust_* for radiative rescaling of dust
470if (dustscaling_mode==2) then
471 if (startphy_file) then
472   call get_field("dust_rad_adjust_prev",dust_rad_adjust_prev,found,indextime)
473   if (.not.found) then
474     write(*,*) "phyetat0: <dust_rad_adjust_prev> not in file; set to 1"
475     dust_rad_adjust_prev(:) = 1
476   endif
477   call get_field("dust_rad_adjust_next",dust_rad_adjust_next,found,indextime)
478   if (.not.found) then
479     write(*,*) "phyetat0: <dust_rad_adjust_next> not in file; set to 1"
480     dust_rad_adjust_next(:) = 1
481   endif
482 else
483   dust_rad_adjust_prev(:)= 0
484   dust_rad_adjust_next(:)= 0
485 endif ! if (startphy_file)
486 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_prev> range:", &
487            minval(dust_rad_adjust_prev), maxval(dust_rad_adjust_prev)
488 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_next> range:", &
489            minval(dust_rad_adjust_next), maxval(dust_rad_adjust_next)
490endif ! of if (dustscaling_mode==2)
491
492! dtau: opacity difference between GCM and dust scenario
493if (startphy_file) then
494   call get_field("dtau",dtau,found,indextime)
495   if (.not.found) then
496     write(*,*) "phyetat0: <dtau> not in file; set to zero"
497     dtau(:) = 0
498   endif
499else
500   dtau(:)= 0
501endif ! if (startphy_file)
502write(*,*) "phyetat0: opacity diff wrt scenario <dtau> range:", &
503            minval(dtau), maxval(dtau)
504
505
506! Sub-grid cloud fraction
507if (startphy_file) then
508   call get_field("totcloudfrac",totcloudfrac,found,indextime)
509   if (.not.found) then
510     write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
511     totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
512   endif
513else
514   totcloudfrac(:)=1.0
515endif ! if (startphy_file)
516write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
517            minval(totcloudfrac), maxval(totcloudfrac)
518
519
520! Max vertical velocity in thermals
521if (startphy_file) then
522   call get_field("wstar",wstar,found,indextime)
523   if (.not.found) then
524     write(*,*) "phyetat0: <wstar> not in file! Set to zero"
525     wstar(:)=0
526   endif
527else
528   wstar(:)=0
529endif ! if (startphy_file)
530write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
531            minval(wstar),maxval(wstar)
532
533
534! Surface temperature :
535if (startphy_file) then !tsurf
536   call get_field("tsurf",tsurf,found,indextime)
537   if (.not.found) then
538     call abort_physic(modname, &
539                "phyetat0: Failed loading <tsurf>",1)
540   endif
541else
542  tsurf(:)=0. ! will be updated afterwards in physiq !
543endif ! of if (startphy_file)
544write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
545            minval(tsurf), maxval(tsurf)
546
547! Surface albedo
548if (startphy_file) then
549   call get_field("albedo",albedo(:,1),found,indextime)
550   if (.not.found) then
551     write(*,*) "phyetat0: Failed loading <albedo>"
552     albedo(:,1)=albedodat(:)
553   endif
554else
555   albedo(:,1)=albedodat(:)
556endif ! of if (startphy_file)
557write(*,*) "phyetat0: Surface albedo <albedo> range:", &
558            minval(albedo(:,1)), maxval(albedo(:,1))
559albedo(:,2)=albedo(:,1)
560
561! Surface emissivity
562if (startphy_file) then
563   call get_field("emis",emis,found,indextime)
564   if (.not.found) then
565     call abort_physic(modname, &
566                "phyetat0: Failed loading <emis>",1)
567   endif
568else
569  ! If no startfi file, use parameter surfemis in def file
570  surfemis=0.95
571  call getin_p("surfemis_without_startfi",surfemis)
572  print*,"surfemis_without_startfi",surfemis
573  emis(:)=surfemis
574endif ! of if (startphy_file)
575write(*,*) "phyetat0: Surface emissivity <emis> range:", &
576            minval(emis), maxval(emis)
577
578
579! surface roughness length (NB: z0 is a common in surfdat_h)
580if (startphy_file) then
581   call get_field("z0",z0,found)
582   if (.not.found) then
583     write(*,*) "phyetat0: Failed loading <z0>"
584     write(*,*) 'will use constant value of z0_default:',z0_default
585     z0(:)=z0_default
586   endif
587else
588   z0(:)=z0_default
589endif ! of if (startphy_file)
590write(*,*) "phyetat0: Surface roughness <z0> range:", &
591            minval(z0), maxval(z0)
592
593
594! pbl wind variance
595if (startphy_file) then
596   call get_field("q2",q2,found,indextime)
597   if (.not.found) then
598     call abort_physic(modname, &
599                "phyetat0: Failed loading <q2>",1)
600   endif
601else
602  q2(:,:)=0.
603endif ! of if (startphy_file)
604write(*,*) "phyetat0: PBL wind variance <q2> range:", &
605            minval(q2), maxval(q2)
606
607! Non-orographic gravity waves
608if (startphy_file) then
609   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
610   if (.not.found) then
611      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
612      du_nonoro_gwd(:,:)=0.
613   endif
614else
615du_nonoro_gwd(:,:)=0.
616endif ! if (startphy_file)
617write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
618write(*,*) " <du_nonoro_gwd> range:", &
619             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
620
621if (startphy_file) then
622   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
623   if (.not.found) then
624      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
625      dv_nonoro_gwd(:,:)=0.
626   endif
627else ! ! if (startphy_file)
628dv_nonoro_gwd(:,:)=0.
629endif ! if (startphy_file)
630write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
631write(*,*) " <dv_nonoro_gwd> range:", &
632             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
633
634! tracer on surface
635if (nq.ge.1) then
636  do iq=1,nq
637    txt=noms(iq)
638    if (txt.eq."h2o_vap") then
639      ! There is no surface tracer for h2o_vap;
640      ! "h2o_ice" should be loaded instead
641      txt="h2o_ice"
642      write(*,*) 'phyetat0: loading surface tracer', &
643                           ' h2o_ice instead of h2o_vap'
644      write(*,*) 'iq = ', iq
645    endif
646
647    if (hdo) then
648    if (txt.eq."hdo_vap") then
649      txt="hdo_ice"
650      write(*,*) 'phyetat0: loading surface tracer', &
651                           ' hdo_ice instead of hdo_vap'
652    endif
653    endif !hdo
654
655    if (startphy_file) then
656      if (txt.eq."co2") then
657        ! We first check if co2ice exist in the startfi file (old way)
658        ! CO2 ice cover
659        call get_field("co2ice",qsurf(:,iq),found,indextime)
660        ! If not, we load the corresponding tracer in qsurf
661        if (.not.found) then
662          call get_field(txt,qsurf(:,iq),found,indextime)
663          if (.not.found) then
664            call abort_physic(modname, &
665                "phyetat0: Failed loading co2ice. there is neither the variable co2ice nor qsurf",1)
666          endif
667        endif
668      else ! (not the co2 tracer)
669        call get_field(txt,qsurf(:,iq),found,indextime)
670        if (.not.found) then
671          write(*,*) "phyetat0: Failed loading <",trim(txt),">"
672          write(*,*) "         ",trim(txt)," is set to zero"
673          qsurf(:,iq)=0.
674        endif
675      endif !endif co2
676    else !(not startphy_file)
677      qsurf(:,iq)=0.
678    endif ! of if (startphy_file)
679    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
680                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
681  enddo ! of do iq=1,nq
682
683    if (txt.eq."co2") then
684      ! We first check if co2ice exist in the startfi file (old way)
685      ! CO2 ice cover
686      if (startphy_file) then
687        call get_field("co2ice",qsurf(:,iq),found,indextime)
688      ! If not, we load the corresponding tracer in qsurf
689        if (.not.found) then
690          call get_field(txt,qsurf(:,iq),found,indextime)
691          if (.not.found) then
692            call abort_physic(modname, &
693                "phyetat0: Failed loading co2ice",1)
694          endif
695        endif
696      else
697       ! If we run without startfile, co2ice is set to 0
698        qsurf(:,iq)=0.
699      endif !if (startphy_file)
700        write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
701            minval(qsurf(:,iq)), maxval(qsurf(:,iq))
702    endif
703
704endif ! of if (nq.ge.1)
705
706if (startphy_file) then
707   call get_field("watercap",watercap,found,indextime)
708   if (.not.found) then
709     write(*,*) "phyetat0: Failed loading <watercap> : ", &
710                          "<watercap> is set to zero"
711     watercap(:)=0
712     write(*,*) 'Now transfer negative surface water ice to', &
713                ' watercap !'
714     if (nq.ge.1) then
715      do iq=1,nq
716       txt=noms(iq)
717       if (txt.eq."h2o_ice") then
718         do ig=1,ngrid
719          if (qsurf(ig,iq).lt.0.0) then
720             watercap(ig) = qsurf(ig,iq)
721             qsurf(ig,iq) = 0.0
722          end if
723         end do
724       endif
725
726       if (txt.eq."hdo_ice") then
727         do ig=1,ngrid
728          if (qsurf(ig,iq).lt.0.0) then
729             qsurf(ig,iq) = 0.0
730          end if
731         end do
732       endif
733
734      enddo
735     endif ! of if (nq.ge.1)
736   endif ! of if (.not.found)
737else
738   watercap(:)=0
739endif ! of if (startphy_file)
740write(*,*) "phyetat0: Surface water ice <watercap> range:", &
741            minval(watercap), maxval(watercap)
742
743if (startphy_file) then
744  ! Call to soil_settings, in order to read soil temperatures,
745  ! as well as thermal inertia and volumetric heat capacity
746  call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
747else
748    flux_geo(:) = 0.
749endif ! of if (startphy_file)
750
751if (startphy_file) then
752   call get_field("watercaptag",watercaptag_tmp,found,indextime)
753   if (.not.found) then
754     write(*,*) "phyetat0: Failed loading <watercaptag> : ", &
755                          "<watercaptag> is set as defined by icelocationmode in surfini.F"
756     watercaptag(:)=.false.
757   else
758     do ig=1,ngrid
759       if(watercaptag_tmp(ig).lt.0.5) then
760          watercaptag(ig)=.false.
761       else
762          watercaptag(ig)=.true.
763       endif
764     enddo
765   endif
766endif !startphy_file
767
768!
769! close file:
770!
771if (startphy_file) call close_startphy
772
773end subroutine phyetat0
774
775
776subroutine ini_tab_controle_dyn_xios(idayref)
777
778      USE comcstfi_h, only: g, mugaz, omeg, rad, rcp
779      USE time_phylmdz_mod, ONLY: hour_ini, daysec, dtphys
780      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
781      IMPLICIT NONE
782
783
784      INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
785
786
787      INTEGER length,l
788      parameter (length = 100)
789      REAL tab_cntrl(length) ! run parameters are stored in this array
790
791      DO l=1,length
792         tab_cntrl(l)=0.
793      ENDDO
794      tab_cntrl(1)  = real(nbp_lon)
795      tab_cntrl(2)  = real(nbp_lat-1)
796      tab_cntrl(3)  = real(nbp_lev)
797      tab_cntrl(4)  = real(idayref)
798      tab_cntrl(5)  = rad
799      tab_cntrl(6)  = omeg
800      tab_cntrl(7)  = g
801      tab_cntrl(8)  = mugaz
802      tab_cntrl(9)  = rcp
803      tab_cntrl(10) = daysec
804      tab_cntrl(11) = dtphys
805
806      tab_cntrl(27) = hour_ini
807
808      tab_cntrl_mod = tab_cntrl
809
810end subroutine ini_tab_controle_dyn_xios
811
812end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.