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

Last change on this file since 2884 was 2884, checked in by romain.vande, 22 months ago

Mars PCM:
Watercaptag is now outputed in the starfi.nc.
There is a new mode : icelocationmode=5 to read watercaptag from startfi and to do as iceloctaionmode=4 if the variable is not present.
RV

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