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

Last change on this file since 2947 was 2919, checked in by llange, 22 months ago

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