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

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

Mars PEM:
New Boolean options for following orbital parameters of ob_ex_lsp.asc: var_obl, var_ex, var_lsp.
If using evol_orbit_pem=true, you can specify which parameter to follow.
True by default: Do you want to change the parameter XXX during the PEM run as prescribed in ob_ex_lsp.asc.
If false, it is set to constant (to the value of the tab_cntrl in the start)
RV
File size: 24.2 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 (startphy_file) then
137 call get_var("def_slope",def_slope,found)
138   if(.not.found) then
139     startphy_slope=.false.
140     write(*,*)'slope_settings: Problem while reading <def_slope>'
141     write(*,*)'Checking that nslope=1'
142     if(nslope.eq.1) then
143      write(*,*)'We take default def_slope and subslope dist'
144      allocate(default_def_slope(nslope+1))
145      default_def_slope(1) = 0.
146      default_def_slope(2) = 0.
147      subslope_dist(:,:)=1.
148     else
149       write(*,*)'Variable def_slope is not present in the start and'
150       write(*,*)'you are trying to run with nslope!=1'
151       write(*,*)'This is not possible, you have to go through newstart'
152     endif
153   else
154     startphy_slope=.true.
155     call get_field("subslope_dist",subslope_dist,found,indextime)
156     if(.not.found) then
157       write(*,*)'slope_settings: Problem while reading <subslope_dist>'
158       write(*,*)'We have to abort.'
159       write(*,*)'Please check that nslope is coherent, you can modify this parameter with newstart'
160       call abort_physic(modname, &
161                "phyetat0: Failed loading <subslope_dist>",1)
162     endif
163    endif
164else ! (no startphy_file)
165     if(nslope.eq.1) then
166      allocate(default_def_slope(nslope+1))
167      default_def_slope(1) = 0.
168      default_def_slope(2) = 0.
169      subslope_dist(:,:)=1.
170    else
171     write(*,*)'Without starfi, lets first run with nslope=1'
172       call abort_physic(modname, &
173                "phyetat0: No startfi and nslope!=1",1)
174   endif
175endif   
176
177do islope=1,nslope
178    def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
179enddo
180 
181DO ig = 1,ngrid
182  sum_dist = 0.
183  DO islope = 1,nslope
184     sum_dist = sum_dist + subslope_dist(ig,islope)
185  ENDDO
186  DO islope = 1,nslope
187    subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
188  ENDDO
189ENDDO
190
191!Now determine the major subslope, ie. the maximal distribution
192
193DO ig=1,ngrid
194  major_slope(ig)=1
195  current_max=subslope_dist(ig,1)
196  DO islope=2,nslope
197    if(subslope_dist(ig,islope).GT.current_max) then
198      major_slope(ig)=islope
199      current_max=subslope_dist(ig,islope)
200    ENDIF
201  ENDDO
202ENDDO
203
204if (startphy_file) then
205   ! Load surface geopotential:
206   call get_field("phisfi",phisfi,found)
207   if (.not.found) then
208     call abort_physic(modname, &
209                "phyetat0: Failed loading <phisfi>",1)
210   endif
211else
212  phisfi(:)=0.
213endif ! of if (startphy_file)
214write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
215               minval(phisfi), maxval(phisfi)
216
217
218if (startphy_file) then
219   ! Load bare ground albedo:
220   call get_field("albedodat",albedodat,found)
221   if (.not.found) then
222     call abort_physic(modname, &
223                "phyetat0: Failed loading <albedodat>",1)
224   endif
225else ! If no startfi file, use parameter surfalbedo in def file
226  surfalbedo=0.1
227  call getin_p("surfalbedo_without_startfi",surfalbedo)
228  print*,"surfalbedo_without_startfi",surfalbedo
229  albedodat(:)=surfalbedo
230endif ! of if (startphy_file)
231write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
232             minval(albedodat), maxval(albedodat)
233
234! ZMEA
235if (startphy_file) then
236   call get_field("ZMEA",zmea,found)
237   if (.not.found) then
238     call abort_physic(modname, &
239                "phyetat0: Failed loading <ZMEA>",1)
240   endif
241else
242  zmea(:)=0.
243endif ! of if (startphy_file)
244write(*,*) "phyetat0: <ZMEA> range:", &
245            minval(zmea), maxval(zmea)
246
247
248! ZSTD
249if (startphy_file) then
250   call get_field("ZSTD",zstd,found)
251   if (.not.found) then
252     call abort_physic(modname, &
253                "phyetat0: Failed loading <ZSTD>",1)
254   endif
255else
256  zstd(:)=0.
257endif ! of if (startphy_file)
258write(*,*) "phyetat0: <ZSTD> range:", &
259            minval(zstd), maxval(zstd)
260
261
262! ZSIG
263if (startphy_file) then
264   call get_field("ZSIG",zsig,found)
265   if (.not.found) then
266     call abort_physic(modname, &
267                "phyetat0: Failed loading <ZSIG>",1)
268   endif
269else
270  zsig(:)=0.
271endif ! of if (startphy_file)
272write(*,*) "phyetat0: <ZSIG> range:", &
273            minval(zsig), maxval(zsig)
274
275
276! ZGAM
277if (startphy_file) then
278   call get_field("ZGAM",zgam,found)
279   if (.not.found) then
280     call abort_physic(modname, &
281                "phyetat0: Failed loading <ZGAM>",1)
282   endif
283else
284  zgam(:)=0.
285endif ! of if (startphy_file)
286write(*,*) "phyetat0: <ZGAM> range:", &
287            minval(zgam), maxval(zgam)
288
289
290! ZTHE
291if (startphy_file) then
292   call get_field("ZTHE",zthe,found)
293   if (.not.found) then
294     call abort_physic(modname, &
295                "phyetat0: Failed loading <ZTHE>",1)
296   endif
297else
298  zthe(:)=0.
299endif ! of if (startphy_file)
300write(*,*) "phyetat0: <ZTHE> range:", &
301             minval(zthe), maxval(zthe)
302
303! hmons
304if (startphy_file) then
305   call get_field("hmons",hmons,found)
306   if (.not.found) then
307     write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
308     if (rdstorm) then
309     call abort_physic(modname, &
310                "phyetat0: Failed loading <hmons>",1)
311     else
312       write(*,*) "will continue anyway..."
313       write(*,*) "because you may not need it."
314       hmons(:)=0.
315     end if ! if (rdstorm)
316   else
317     do ig=1,ngrid
318       if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
319     enddo
320   endif ! (.not.found)
321else
322   hmons(:)=0.
323endif ! if (startphy_file)
324write(*,*) "phyetat0: <hmons> range:", &
325            minval(hmons), maxval(hmons)
326
327
328! summit
329if (startphy_file) then
330   call get_field("summit",summit,found)
331   if (.not.found) then
332     write(*,*) "WARNING: phyetat0: Failed loading <summit>"
333     if (rdstorm) then
334     call abort_physic(modname, &
335                "phyetat0: Failed loading <summit>",1)
336     else
337       write(*,*) "will continue anyway..."
338       write(*,*) "because you may not need it."
339       summit(:)=0.
340     end if
341   else
342     do ig=1,ngrid
343       if (summit(ig) .eq. -999999.)  summit(ig)=0.
344     enddo
345   endif ! if (.not.found)
346else
347   summit(:)=0. 
348endif ! if (startphy_file)
349write(*,*) "phyetat0: <summit> range:", &
350            minval(summit), maxval(summit)
351
352
353! base
354if (startphy_file) then
355   call get_field("base",base,found)
356   if (.not.found) then
357     write(*,*) "WARNING: phyetat0: Failed loading <base>"
358     if (rdstorm) then
359     call abort_physic(modname, &
360                "phyetat0: Failed loading <base>",1)
361     else
362       write(*,*) "will continue anyway..."
363       write(*,*) "because you may not need it."
364       base(:)=0.
365     end if
366   else
367     do ig=1,ngrid
368       if (base(ig) .eq. -999999.)  base(ig)=0.
369     enddo
370   endif ! if(.not.found)
371else
372   base(:)=0.
373endif ! if (startphy_file)
374write(*,*) "phyetat0: <base> range:", &
375            minval(base), maxval(base)
376
377! Time axis
378! obtain timestart from run.def
379timestart=-9999 ! default value
380call getin_p("timestart",timestart)
381if (startphy_file) then
382   found=inquire_dimension("Time")
383   if (.not.found) then
384     indextime = 1
385     write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
386   else
387     write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
388     timelen=inquire_dimension_length("Time")
389     allocate(time(timelen))
390     ! load "Time" array:
391     call get_var("Time",time,found)
392     if (.not.found) then
393     call abort_physic(modname, &
394                "phyetat0: Failed loading <Time>",1)
395     endif
396     ! seclect the desired time index
397     IF (timestart .lt. 0) THEN  ! default: we use the last time value
398       indextime = timelen
399     ELSE  ! else we look for the desired value in the time axis
400       indextime = 0
401       DO i=1,timelen
402         IF (abs(time(i) - timestart) .lt. 0.01) THEN
403           indextime = i
404           EXIT
405         ENDIF
406       ENDDO
407       IF (indextime .eq. 0) THEN
408         PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
409         PRINT*, "Stored times are:"
410         DO i=1,timelen
411            PRINT*, time(i)
412         ENDDO
413         call abort_physic(modname,"phyetat0: Time error",1)
414       ENDIF
415     ENDIF ! of IF (timestart .lt. 0)
416     ! In startfi the absolute date is day_ini + time0 + time
417     ! For now on, in the GCM physics, it is day_ini + time0
418     time0 = time(indextime) + time0
419     day_ini = day_ini + INT(time0)
420     time0 = time0 - INT(time0)   
421     PRINT*, "phyetat0: Selected time ",time(indextime), &
422             " at index ",indextime
423     DEALLOCATE(time)
424   endif ! of if Time not found in file
425
426   call ini_tab_controle_dyn_xios(day_ini)
427
428else
429  indextime = 1
430endif ! if (startphy_file)
431
432! Dust conversion factor
433if (startphy_file) then
434   call get_field("tauscaling",tauscaling,found,indextime)
435   if (.not.found) then
436     write(*,*) "phyetat0: <tauscaling> not in file"
437     tauscaling(:) = 1
438   endif
439else
440   tauscaling(:) = 1
441endif ! if (startphy_file)
442write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
443            minval(tauscaling), maxval(tauscaling)
444
445! dust_rad_adjust_* for radiative rescaling of dust
446if (dustscaling_mode==2) then
447 if (startphy_file) then
448   call get_field("dust_rad_adjust_prev",dust_rad_adjust_prev,found,indextime)
449   if (.not.found) then
450     write(*,*) "phyetat0: <dust_rad_adjust_prev> not in file; set to 1"
451     dust_rad_adjust_prev(:) = 1
452   endif
453   call get_field("dust_rad_adjust_next",dust_rad_adjust_next,found,indextime)
454   if (.not.found) then
455     write(*,*) "phyetat0: <dust_rad_adjust_next> not in file; set to 1"
456     dust_rad_adjust_next(:) = 1
457   endif
458 else
459   dust_rad_adjust_prev(:)= 0
460   dust_rad_adjust_next(:)= 0
461 endif ! if (startphy_file)
462 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_prev> range:", &
463            minval(dust_rad_adjust_prev), maxval(dust_rad_adjust_prev)
464 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_next> range:", &
465            minval(dust_rad_adjust_next), maxval(dust_rad_adjust_next)
466endif ! of if (dustscaling_mode==2)
467
468! dtau: opacity difference between GCM and dust scenario
469if (startphy_file) then
470   call get_field("dtau",dtau,found,indextime)
471   if (.not.found) then
472     write(*,*) "phyetat0: <dtau> not in file; set to zero"
473     dtau(:) = 0
474   endif
475else
476   dtau(:)= 0
477endif ! if (startphy_file)
478write(*,*) "phyetat0: opacity diff wrt scenario <dtau> range:", &
479            minval(dtau), maxval(dtau)
480
481
482! Sub-grid cloud fraction
483if (startphy_file) then
484   call get_field("totcloudfrac",totcloudfrac,found,indextime)
485   if (.not.found) then
486     write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
487     totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
488   endif
489else
490   totcloudfrac(:)=1.0
491endif ! if (startphy_file)
492write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
493            minval(totcloudfrac), maxval(totcloudfrac)
494
495
496! Max vertical velocity in thermals
497if (startphy_file) then
498   call get_field("wstar",wstar,found,indextime)
499   if (.not.found) then
500     write(*,*) "phyetat0: <wstar> not in file! Set to zero"
501     wstar(:)=0
502   endif
503else
504   wstar(:)=0
505endif ! if (startphy_file)
506write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
507            minval(wstar),maxval(wstar)
508
509
510! Surface temperature :
511if (startphy_file) then !tsurf
512   call get_field("tsurf",tsurf,found,indextime)
513   if (.not.found) then
514     call abort_physic(modname, &
515                "phyetat0: Failed loading <tsurf>",1)
516   endif
517else
518  tsurf(:)=0. ! will be updated afterwards in physiq !
519endif ! of if (startphy_file)
520write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
521            minval(tsurf), maxval(tsurf)
522
523! Surface albedo
524if (startphy_file) then
525   call get_field("albedo",albedo(:,1),found,indextime)
526   if (.not.found) then
527     write(*,*) "phyetat0: Failed loading <albedo>"
528     albedo(:,1)=albedodat(:)
529   endif
530else
531   albedo(:,1)=albedodat(:)
532endif ! of if (startphy_file)
533write(*,*) "phyetat0: Surface albedo <albedo> range:", &
534            minval(albedo(:,1)), maxval(albedo(:,1))
535albedo(:,2)=albedo(:,1)
536
537! Surface emissivity
538if (startphy_file) then
539   call get_field("emis",emis,found,indextime)
540   if (.not.found) then
541     call abort_physic(modname, &
542                "phyetat0: Failed loading <emis>",1)
543   endif
544else
545  ! If no startfi file, use parameter surfemis in def file
546  surfemis=0.95
547  call getin_p("surfemis_without_startfi",surfemis)
548  print*,"surfemis_without_startfi",surfemis
549  emis(:)=surfemis
550endif ! of if (startphy_file)
551write(*,*) "phyetat0: Surface emissivity <emis> range:", &
552            minval(emis), maxval(emis)
553
554
555! surface roughness length (NB: z0 is a common in surfdat_h)
556if (startphy_file) then
557   call get_field("z0",z0,found)
558   if (.not.found) then
559     write(*,*) "phyetat0: Failed loading <z0>"
560     write(*,*) 'will use constant value of z0_default:',z0_default
561     z0(:)=z0_default
562   endif
563else
564   z0(:)=z0_default
565endif ! of if (startphy_file)
566write(*,*) "phyetat0: Surface roughness <z0> range:", &
567            minval(z0), maxval(z0)
568
569
570! pbl wind variance
571if (startphy_file) then
572   call get_field("q2",q2,found,indextime)
573   if (.not.found) then
574     call abort_physic(modname, &
575                "phyetat0: Failed loading <q2>",1)
576   endif
577else
578  q2(:,:)=0.
579endif ! of if (startphy_file)
580write(*,*) "phyetat0: PBL wind variance <q2> range:", &
581            minval(q2), maxval(q2)
582
583! Non-orographic gravity waves
584if (startphy_file) then
585   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
586   if (.not.found) then
587      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
588      du_nonoro_gwd(:,:)=0.
589   endif
590else
591du_nonoro_gwd(:,:)=0.
592endif ! if (startphy_file)
593write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
594write(*,*) " <du_nonoro_gwd> range:", &
595             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
596
597if (startphy_file) then
598   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
599   if (.not.found) then
600      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
601      dv_nonoro_gwd(:,:)=0.
602   endif
603else ! ! if (startphy_file)
604dv_nonoro_gwd(:,:)=0.
605endif ! if (startphy_file)
606write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
607write(*,*) " <dv_nonoro_gwd> range:", &
608             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
609
610! tracer on surface
611if (nq.ge.1) then
612  do iq=1,nq
613    txt=noms(iq)
614    if (txt.eq."h2o_vap") then
615      ! There is no surface tracer for h2o_vap;
616      ! "h2o_ice" should be loaded instead
617      txt="h2o_ice"
618      write(*,*) 'phyetat0: loading surface tracer', &
619                           ' h2o_ice instead of h2o_vap'
620      write(*,*) 'iq = ', iq
621    endif
622
623    if (hdo) then
624    if (txt.eq."hdo_vap") then
625      txt="hdo_ice"
626      write(*,*) 'phyetat0: loading surface tracer', &
627                           ' hdo_ice instead of hdo_vap'
628    endif
629    endif !hdo
630
631    if (startphy_file) then
632      if (txt.eq."co2") then
633        ! We first check if co2ice exist in the startfi file (old way)
634        ! CO2 ice cover
635        call get_field("co2ice",qsurf(:,iq),found,indextime)
636        ! If not, we load the corresponding tracer in qsurf
637        if (.not.found) then
638          call get_field(txt,qsurf(:,iq),found,indextime)
639          if (.not.found) then
640            call abort_physic(modname, &
641                "phyetat0: Failed loading co2ice. there is neither the variable co2ice nor qsurf",1)
642          endif
643        endif
644      else ! (not the co2 tracer)
645        call get_field(txt,qsurf(:,iq),found,indextime)
646        if (.not.found) then
647          write(*,*) "phyetat0: Failed loading <",trim(txt),">"
648          write(*,*) "         ",trim(txt)," is set to zero"
649          qsurf(:,iq)=0.
650        endif
651      endif !endif co2
652    else !(not startphy_file)
653      qsurf(:,iq)=0.
654    endif ! of if (startphy_file)
655    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
656                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
657  enddo ! of do iq=1,nq
658
659    if (txt.eq."co2") then
660      ! We first check if co2ice exist in the startfi file (old way)
661      ! CO2 ice cover
662      if (startphy_file) then
663        call get_field("co2ice",qsurf(:,iq),found,indextime)
664      ! If not, we load the corresponding tracer in qsurf
665        if (.not.found) then
666          call get_field(txt,qsurf(:,iq),found,indextime)
667          if (.not.found) then
668            call abort_physic(modname, &
669                "phyetat0: Failed loading co2ice",1)
670          endif
671        endif
672      else
673       ! If we run without startfile, co2ice is set to 0
674        qsurf(:,iq)=0.
675      endif !if (startphy_file)
676        write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
677            minval(qsurf(:,iq)), maxval(qsurf(:,iq))
678    endif
679
680endif ! of if (nq.ge.1)
681
682if (startphy_file) then
683   call get_field("watercap",watercap,found,indextime)
684   if (.not.found) then
685     write(*,*) "phyetat0: Failed loading <watercap> : ", &
686                          "<watercap> is set to zero"
687     watercap(:)=0
688     write(*,*) 'Now transfer negative surface water ice to', &
689                ' watercap !'
690     if (nq.ge.1) then
691      do iq=1,nq
692       txt=noms(iq)
693       if (txt.eq."h2o_ice") then
694         do ig=1,ngrid
695          if (qsurf(ig,iq).lt.0.0) then
696             watercap(ig) = qsurf(ig,iq)
697             qsurf(ig,iq) = 0.0
698          end if
699         end do
700       endif
701
702       if (txt.eq."hdo_ice") then
703         do ig=1,ngrid
704          if (qsurf(ig,iq).lt.0.0) then
705             qsurf(ig,iq) = 0.0
706          end if
707         end do
708       endif
709
710      enddo
711     endif ! of if (nq.ge.1)
712   endif ! of if (.not.found)
713else
714   watercap(:)=0
715endif ! of if (startphy_file)
716write(*,*) "phyetat0: Surface water ice <watercap> range:", &
717            minval(watercap), maxval(watercap)
718
719if (startphy_file) then
720  ! Call to soil_settings, in order to read soil temperatures,
721  ! as well as thermal inertia and volumetric heat capacity
722  call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
723else
724    flux_geo(:) = 0.
725endif ! of if (startphy_file)
726
727if (startphy_file) then
728   call get_field("watercaptag",watercaptag_tmp,found,indextime)
729   if (.not.found) then
730     write(*,*) "phyetat0: Failed loading <watercaptag> : ", &
731                          "<watercaptag> is set as defined by icelocationmode in surfini.F"
732     watercaptag(:)=.false.
733   else
734     do ig=1,ngrid
735       if(watercaptag_tmp(ig).lt.0.5) then
736          watercaptag(ig)=.false.
737       else
738          watercaptag(ig)=.true.
739       endif
740     enddo
741   endif
742endif !startphy_file
743
744!
745! close file:
746!
747if (startphy_file) call close_startphy
748
749end subroutine phyetat0
750
751
752subroutine ini_tab_controle_dyn_xios(idayref)
753
754      USE comcstfi_h, only: g, mugaz, omeg, rad, rcp
755      USE time_phylmdz_mod, ONLY: hour_ini, daysec, dtphys
756      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
757      IMPLICIT NONE
758
759
760      INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
761
762
763      INTEGER length,l
764      parameter (length = 100)
765      REAL tab_cntrl(length) ! run parameters are stored in this array
766
767      DO l=1,length
768         tab_cntrl(l)=0.
769      ENDDO
770      tab_cntrl(1)  = real(nbp_lon)
771      tab_cntrl(2)  = real(nbp_lat-1)
772      tab_cntrl(3)  = real(nbp_lev)
773      tab_cntrl(4)  = real(idayref)
774      tab_cntrl(5)  = rad
775      tab_cntrl(6)  = omeg
776      tab_cntrl(7)  = g
777      tab_cntrl(8)  = mugaz
778      tab_cntrl(9)  = rcp
779      tab_cntrl(10) = daysec
780      tab_cntrl(11) = dtphys
781
782      tab_cntrl(27) = hour_ini
783
784      tab_cntrl_mod = tab_cntrl
785
786end subroutine ini_tab_controle_dyn_xios
787
788end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.