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

Last change on this file since 2997 was 2994, checked in by llange, 18 months ago

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