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

Last change on this file since 3906 was 3903, checked in by jbclement, 6 months ago

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