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

Last change on this file since 3064 was 3037, checked in by jbclement, 22 months ago

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