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

Last change on this file since 2443 was 2417, checked in by emillour, 5 years ago

Mars GCM:
Add a new scheme to handle the dust and its radiative impact. Triggered using
a new flag dustscaling_mode=2 (dustscaling_mod=0: no rescaling at all, and
dustscaling_mode=1: full rescaling using tauscaling, GCMv5.3 style). Rescaling
is then only done on the radiative impact (see dust_scaling_mod.F90) of dust.
Moreover the scaling factor "dust_rad_adjust" is evaluated using the target dust
scenario opacity for the next sol and left to evolve linearly until then to not
impose the diurnal evolution of dust.
In practice, main changes or additions in the code are:

  • renamed flag "tauscaling_mode" as "dustscaling_mode"
  • moved parameter "t_scenario_sol" to "dust_param_mod"
  • adapted "compute_dustscaling" routine in "dust_scaling_mod"
  • added module "dust_rad_adjust_mod"
  • 2D fields "dust_rad_adjust_prev" and "dust_rad_adjust_next" required to compute coefficient "dust_rad_adjust" need to be stored in (re)startfi files

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