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

Last change on this file since 3026 was 2999, checked in by llange, 18 months ago

Mars PCM
Include perenial_co2ice (equivalent of watercap) to distinguich between CO2 frost and perenial CO2 ice for paleoclimate studies.
When no frost is present and we dig into perenial ice, the surface albedo is changed. The albedo for seasonal ice is set to 0.65, and the perenial ice albedo can be fixed in the callphys.def. I recommand values between 0.8 and 0.9.
To use this, paleoclimate must be set to True and TESalbedo to false in the callphys.def. Else, it runs as usual with TES albedo
LL
File size: 26.7 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.
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
670
671    if (txt.eq."co2") then
672      ! We first check if co2ice exist in the startfi file (old way)
673      ! CO2 ice cover
674      if (startphy_file) then
675        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
676      ! If not, we load the corresponding tracer in qsurf
677        if (.not.found) then
678          call get_field(txt,qsurf(:,iq,:),found,indextime)
679          if (.not.found) then
680            call abort_physic(modname, &
681                "phyetat0: Failed loading co2ice",1)
682          endif
683        endif
684      else
685       ! If we run without startfile, co2ice is set to 0
686        qsurf(:,iq,:)=0.
687      endif !if (startphy_file)
688        write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
689            minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
690    endif
691
692endif ! of if (nq.ge.1)
693
694if (startphy_file) then
695   call get_field("watercap",watercap,found,indextime)
696   if (.not.found) then
697     write(*,*) "phyetat0: Failed loading <watercap> : ", &
698                          "<watercap> is set to zero"
699     watercap(:,:)=0
700     write(*,*) 'Now transfer negative surface water ice to', &
701                ' watercap !'
702     if (nq.ge.1) then
703      do iq=1,nq
704       txt=noms(iq)
705       if (txt.eq."h2o_ice") then
706         do ig=1,ngrid
707          do islope=1,nslope
708            if (qsurf(ig,iq,islope).lt.0.0) then
709              watercap(ig,islope) = qsurf(ig,iq,islope)
710              qsurf(ig,iq,islope) = 0.0
711            end if
712          enddo
713         end do
714       endif
715
716       if (txt.eq."hdo_ice") then
717         do ig=1,ngrid
718          do islope=1,nslope
719            if (qsurf(ig,iq,islope).lt.0.0) then
720              qsurf(ig,iq,islope) = 0.0
721            end if
722          enddo
723         end do
724       endif
725
726      enddo
727     endif ! of if (nq.ge.1)
728   endif ! of if (.not.found)
729else
730   watercap(:,:)=0
731endif ! of if (startphy_file)
732write(*,*) "phyetat0: Surface water ice <watercap> range:", &
733            minval(watercap), maxval(watercap)
734
735if (startphy_file) then
736  ! Call to soil_settings, in order to read soil temperatures,
737  ! as well as thermal inertia and volumetric heat capacity
738  call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
739else
740    flux_geo(:,:) = 0.
741endif ! of if (startphy_file)
742
743if (startphy_file) then
744   call get_field("watercaptag",watercaptag_tmp,found,indextime)
745   if (.not.found) then
746     write(*,*) "phyetat0: Failed loading <watercaptag> : ", &
747                          "<watercaptag> is set as defined by icelocationmode in surfini.F"
748     watercaptag(:)=.false.
749   else
750     do ig=1,ngrid
751       if(watercaptag_tmp(ig).lt.0.5) then
752          watercaptag(ig)=.false.
753       else
754          watercaptag(ig)=.true.
755       endif
756     enddo
757   endif
758endif !startphy_file
759
760
761if (paleoclimate) then
762  do iq=1,nq
763   txt=noms(iq)
764   if (txt.eq."co2") igcm_co2_tmp = iq
765  enddo
766  if (startphy_file) then
767! Depth of H2O lag
768   call get_field("lag_h2o_ice",lag_h2o_ice,found,indextime)
769   if (.not.found) then
770     write(*,*) "phyetat0: Failed loading <lag_h2o_ice> : ", &
771                          "<lag_h2o_ice> is set as -1 (no subsurface ice)"
772     lag_h2o_ice(:,:) = -1.
773   endif
774
775!  Depth of CO2 lag
776   call get_field("lag_co2_ice",lag_co2_ice,found,indextime)
777   if (.not.found) then
778     write(*,*) "phyetat0: Failed loading <lag_co2_ice> : ", &
779                          "<lag_co2_ice> is set as -1 (no subsurface ice)"
780     lag_co2_ice(:,:) = -1.
781   endif
782
783! Perenial CO2 ice   
784   call get_field("perenial_co2ice",perenial_co2ice,found,indextime)
785   if (.not.found) then
786     write(*,*) "phyetat0: Failed loading <perenial_co2ice> : ", &
787                "<perenial_co2ice> is set as 10m at the South Pole"
788     perenial_co2ice(:,:) = 0.
789     if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
790       DO islope = 1,nslope
791        perenial_co2ice(ngrid,islope) = 10*1.6e3 ! 10m which is convert to kg/m^2
792        qsurf(ngrid,igcm_co2_tmp,islope) = qsurf(ngrid-1,igcm_co2_tmp,islope) + & ! perenial ice + frost
793                                        perenial_co2ice(ngrid,islope)
794       ENDDO
795     endif
796   endif ! notfound
797  else ! no startfiphyle
798     lag_h2o_ice(:,:) = -1.
799     lag_co2_ice(:,:) = -1.
800     perenial_co2ice(:,:) = 0.
801     if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
802       DO islope = 1,nslope
803        perenial_co2ice(ngrid,islope) = 10*1.6e3 ! 10m which is convert to kg/m^2
804        qsurf(ngrid,igcm_co2_tmp,islope) = qsurf(ngrid-1,igcm_co2_tmp,islope) + & ! perenial ice + frost
805                                        perenial_co2ice(ngrid,islope)
806       ENDDO
807     endif
808  endif !startphy_file
809else
810  lag_h2o_ice(:,:) = -1.
811  lag_co2_ice(:,:) = -1.
812  perenial_co2ice(:,:) = 0.
813endif !paleoclimate
814
815! close file:
816!
817if (startphy_file) call close_startphy
818
819end subroutine phyetat0
820
821
822subroutine ini_tab_controle_dyn_xios(idayref)
823
824      USE comcstfi_h, only: g, mugaz, omeg, rad, rcp
825      USE time_phylmdz_mod, ONLY: hour_ini, daysec, dtphys
826      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
827      IMPLICIT NONE
828
829
830      INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
831
832
833      INTEGER length,l
834      parameter (length = 100)
835      REAL tab_cntrl(length) ! run parameters are stored in this array
836
837      DO l=1,length
838         tab_cntrl(l)=0.
839      ENDDO
840      tab_cntrl(1)  = real(nbp_lon)
841      tab_cntrl(2)  = real(nbp_lat-1)
842      tab_cntrl(3)  = real(nbp_lev)
843      tab_cntrl(4)  = real(idayref)
844      tab_cntrl(5)  = rad
845      tab_cntrl(6)  = omeg
846      tab_cntrl(7)  = g
847      tab_cntrl(8)  = mugaz
848      tab_cntrl(9)  = rcp
849      tab_cntrl(10) = daysec
850      tab_cntrl(11) = dtphys
851
852      tab_cntrl(27) = hour_ini
853
854      tab_cntrl_mod = tab_cntrl
855
856end subroutine ini_tab_controle_dyn_xios
857
858end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.