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

Last change on this file since 4022 was 3974, checked in by jbclement, 3 months ago

Mars PCM:
Simplification of the orbital parameters initialization (resolve ticket #109): removing redundant actions in subroutine 'iniorbit' in regard of what is done in "tabfi.F" and moving the subroutine 'lsp2solp' from "tabfi.F" to "planete_h.F90" + some cleanings.
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, CLFvarying, CLFfixval
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 water ice cloud fraction
486totcloudfrac(:)=1.0 ! default value; no cloud fraction
487if (startphy_file) then
488  if (CLFvarying) then
489   call get_field("totcloudfrac",totcloudfrac,found,indextime)
490   if (.not.found) then
491     write(*,*) "phyetat0: <totcloudfrac> not in file; set to ",CLFfixval
492   endif
493   totcloudfrac(:)=CLFfixval
494  endif ! of if (CLFvarying)
495endif ! if (startphy_file)
496write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
497            minval(totcloudfrac), maxval(totcloudfrac)
498
499! Max vertical velocity in thermals
500if (startphy_file) then
501   call get_field("wstar",wstar,found,indextime)
502   if (.not.found) then
503     write(*,*) "phyetat0: <wstar> not in file! Set to zero"
504     wstar(:)=0
505   endif
506else
507   wstar(:)=0
508endif ! if (startphy_file)
509write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
510            minval(wstar),maxval(wstar)
511
512! 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! pbl wind variance
576if (startphy_file) then
577   call get_field("q2",q2,found,indextime)
578   if (.not.found) then
579     call abort_physic(modname, &
580                "phyetat0: Failed loading <q2>",1)
581   endif
582else
583  q2(:,:)=0.
584endif ! of if (startphy_file)
585write(*,*) "phyetat0: PBL wind variance <q2> range:", &
586            minval(q2), maxval(q2)
587
588! Non-orographic gravity waves
589if (startphy_file) then
590   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
591   if (.not.found) then
592      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
593      du_nonoro_gwd(:,:)=0.
594   endif
595else
596du_nonoro_gwd(:,:)=0.
597endif ! if (startphy_file)
598write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
599write(*,*) " <du_nonoro_gwd> range:", &
600             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
601
602if (startphy_file) then
603   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
604   if (.not.found) then
605      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
606      dv_nonoro_gwd(:,:)=0.
607   endif
608else ! ! if (startphy_file)
609dv_nonoro_gwd(:,:)=0.
610endif ! if (startphy_file)
611write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
612write(*,*) " <dv_nonoro_gwd> range:", &
613             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
614if (startphy_file) then
615   call get_field("du_eddymix_gwd",du_eddymix_gwd,found,indextime)
616   if (.not.found) then
617      write(*,*) "phyetat0: <du_eddymix_gwd> not in file"
618      du_eddymix_gwd(:,:)=0.
619   endif
620else
621du_eddymix_gwd(:,:)=0.
622endif ! if (startphy_file)
623write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW mixing"
624write(*,*) " <du_eddymix_gwd> range:", &
625             minval(du_eddymix_gwd), maxval(du_eddymix_gwd)
626
627if (startphy_file) then
628   call get_field("dh_eddymix_gwd",dh_eddymix_gwd,found,indextime)
629   if (.not.found) then
630      write(*,*) "phyetat0: <dh_eddymix_gwd> not in file"
631      dh_eddymix_gwd(:,:)=0.
632   endif
633else
634dh_eddymix_gwd(:,:)=0.
635endif ! if (startphy_file)
636write(*,*) "phyetat0: Memory of PT tendency due to non-orographic GW mixing"
637write(*,*) " <dh_eddymix_gwd> range:", &
638             minval(dh_eddymix_gwd), maxval(dh_eddymix_gwd)
639
640
641if (startphy_file) then
642   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
643   if (.not.found) then
644      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
645      dv_nonoro_gwd(:,:)=0.
646   endif
647else ! ! if (startphy_file)
648dv_nonoro_gwd(:,:)=0.
649endif ! if (startphy_file)
650write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
651write(*,*) " <dv_nonoro_gwd> range:", &
652             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
653
654if (startphy_file) then
655   call get_field("dv_eddymix_gwd",dv_eddymix_gwd,found,indextime)
656   if (.not.found) then
657      write(*,*) "phyetat0: <dv_eddymix_gwd> not in file"
658      dv_eddymix_gwd(:,:)=0.
659   endif
660else ! ! if (startphy_file)
661dv_eddymix_gwd(:,:)=0.
662endif ! if (startphy_file)
663write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
664write(*,*) " <dv_eddymix_gwd> range:", &
665             minval(dv_eddymix_gwd), maxval(dv_eddymix_gwd)
666
667if (startphy_file) then
668   call get_field("dq_eddymix_gwd",dq_eddymix_gwd,found,indextime)
669   if (.not.found) then
670      write(*,*) "phyetat0: <dq_eddymix_gwd> not in file"
671      dq_eddymix_gwd(:,:,:)=0.
672   endif
673else ! ! if (startphy_file)
674dq_eddymix_gwd(:,:,:)=0.
675endif ! if (startphy_file)
676write(*,*) "phyetat0: Memory of tracers tendency due to non-orographic GW mixing"
677write(*,*) " <dq_eddymix_gwd> range:", &
678             minval(dq_eddymix_gwd), maxval(dq_eddymix_gwd)
679
680     
681!if (startphy_file) then
682!   call get_field("dr_depflux_gwd",dr_depflux_gwd,found,indextime)
683!   if (.not.found) then
684!      write(*,*) "phyetat0: <dr_depflux_gwd> not in file"
685!      dr_depflux_gwd(:,:,:)=0.
686!   endif
687!else ! ! if (startphy_file)
688!!dr_depflux_gwd(:,:,:)=0.
689!endif ! if (startphy_file)
690!write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
691!write(*,*) " <dr_depflux_gwd> range:", &
692!             minval(dr_depflux_gwd), maxval(dr_depflux_gwd)
693
694if (startphy_file) then
695   call get_field("de_eddymix_rto",de_eddymix_rto,found,indextime)
696   if (.not.found) then
697      write(*,*) "phyetat0: <de_eddymix_rto> not in file"
698      de_eddymix_rto(:,:)=0.
699   endif
700else ! ! if (startphy_file)
701de_eddymix_rto(:,:)=0.
702endif ! if (startphy_file)
703write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
704write(*,*) " <de_eddymix_rto> range:", &
705             minval(de_eddymix_rto), maxval(de_eddymix_rto)
706
707if (startphy_file) then
708   call get_field("df_eddymix_flx ",df_eddymix_flx ,found,indextime)
709   if (.not.found) then
710      write(*,*) "phyetat0: <df_eddymix_flx > not in file"
711      df_eddymix_flx (:,:)=0.
712   endif
713else ! ! if (startphy_file)
714df_eddymix_flx (:,:)=0.
715endif ! if (startphy_file)
716write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
717write(*,*) " <df_eddymix_flx > range:", &
718             minval(df_eddymix_flx ), maxval(df_eddymix_flx )
719
720! tracer on surface
721if (nq.ge.1) then
722  do iq=1,nq
723    txt=noms(iq)
724    if (txt.eq."h2o_vap") then
725      ! There is no surface tracer for h2o_vap;
726      ! "h2o_ice" should be loaded instead
727      txt="h2o_ice"
728      write(*,*) 'phyetat0: loading surface tracer', &
729                           ' h2o_ice instead of h2o_vap'
730      write(*,*) 'iq = ', iq
731    endif
732
733    if (hdo) then
734      if (txt.eq."hdo_vap") then
735        txt="hdo_ice"
736        write(*,*) 'phyetat0: loading surface tracer', &
737                           ' hdo_ice instead of hdo_vap'
738      endif
739    endif !hdo
740
741    if (startphy_file) then
742      if (txt.eq."co2") then
743        ! We first check if co2ice exist in the startfi file (old way)
744        ! CO2 ice cover
745        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
746        ! If not, we load the corresponding tracer in qsurf
747        if (.not.found) then
748          call get_field(txt,qsurf(:,iq,:),found,indextime)
749          if (.not.found) then
750            call abort_physic(modname, &
751                "phyetat0: Failed loading co2ice. There is neither the variable co2ice nor qsurf!",1)
752          endif
753        endif
754      else ! (not the co2 tracer)
755        call get_field(txt,qsurf(:,iq,:),found,indextime)
756        if (.not.found) then
757          write(*,*) "phyetat0: Failed loading <",trim(txt),">"
758          write(*,*) "         ",trim(txt)," is set to zero"
759          qsurf(:,iq,:)=0.
760        endif
761      endif !endif co2
762    else !(not startphy_file)
763      qsurf(:,iq,:)=0. ! co2ice is set to 0
764    endif ! of if (startphy_file)
765    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
766                 minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
767  enddo ! of do iq=1,nq
768endif ! of if (nq.ge.1)
769
770if (startphy_file) then
771   call get_field("watercap",watercap,found,indextime)
772   if (.not.found) then
773     write(*,*) "phyetat0: Failed loading <watercap> : ", &
774                          "<watercap> is set to zero"
775     watercap(:,:)=0
776     write(*,*) 'Now transfer negative surface water ice to', &
777                ' watercap !'
778     if (nq.ge.1) then
779      do iq=1,nq
780       txt=noms(iq)
781       if (txt.eq."h2o_ice") then
782         do ig=1,ngrid
783          do islope=1,nslope
784            if (qsurf(ig,iq,islope).lt.0.0) then
785              watercap(ig,islope) = qsurf(ig,iq,islope)
786              qsurf(ig,iq,islope) = 0.0
787            end if
788          enddo
789         end do
790       endif
791
792       if (txt.eq."hdo_ice") then
793         do ig=1,ngrid
794          do islope=1,nslope
795            if (qsurf(ig,iq,islope).lt.0.0) then
796              qsurf(ig,iq,islope) = 0.0
797            end if
798          enddo
799         end do
800       endif
801
802      enddo
803     endif ! of if (nq.ge.1)
804   endif ! of if (.not.found)
805else
806   watercap(:,:)=0
807endif ! of if (startphy_file)
808write(*,*) "phyetat0: Surface water ice <watercap> range:", &
809            minval(watercap), maxval(watercap)
810
811if (startphy_file) then
812  ! Call to soil_settings, in order to read soil temperatures,
813  ! as well as thermal inertia and volumetric heat capacity
814  call soil_settings(nid_start,ngrid,nsoil,nqsoil,tsurf,tsoil,qsoil,indextime)
815else
816    flux_geo(:,:) = 0.
817endif ! of if (startphy_file)
818
819if (startphy_file) then
820   call get_field("watercaptag",watercaptag_tmp,found,indextime)
821   if (.not.found) then
822     write(*,*) "phyetat0: Failed loading <watercaptag> : ", &
823                          "<watercaptag> is set as defined by icelocationmode in surfini.F"
824     icelocationmode = 4
825   else
826     do ig=1,ngrid
827       if(watercaptag_tmp(ig) < 0.5) then
828          watercaptag(ig) = .false.
829       else
830          watercaptag(ig) = .true.
831       endif
832     enddo
833   endif
834else
835    watercaptag(:) = .false.
836endif ! of if (startphy_file)
837
838if (paleoclimate) then
839  do iq=1,nq
840   txt=noms(iq)
841   if (txt.eq."co2") igcm_co2_tmp = iq
842  enddo
843
844  if (startphy_file) then
845! Depth of H2O lag
846   call get_field("h2o_ice_depth",h2o_ice_depth,found,indextime)
847   if (.not.found) then
848     write(*,*) "phyetat0: Failed loading <h2o_ice_depth> : ", &
849                          "<h2o_ice_depth> is set as -1 (no subsurface ice)"
850     h2o_ice_depth(:,:) = -1.
851   endif
852   
853! Total flux with SSI
854   call get_field("zdqsdif_ssi_tot",zdqsdif_ssi_tot,found,indextime)
855   if (.not.found) then
856     write(*,*) "phyetat0: Failed loading <zdqsdif_ssi_tot> : ", &
857                          "<zdqsdif_ssi_tot> is set as -1 (no subsurface ice)"
858     zdqsdif_ssi_tot(:,:) = 0.
859   endif
860
861! Diffusion coeficent
862   call get_field("d_coef",d_coef,found,indextime)
863   if (.not.found) then
864     write(*,*) "phyetat0: Failed loading <d_coef> : ", &
865                          "<d_coef> is set as 4e-4 (defualt value)"
866     d_coef(:,:) = 4e-4
867   endif
868
869! Depth of CO2 lag
870   call get_field("lag_co2_ice",lag_co2_ice,found,indextime)
871   if (.not.found) then
872     write(*,*) "phyetat0: Failed loading <lag_co2_ice> : ", &
873                          "<lag_co2_ice> is set as -1 (no subsurface ice)"
874     lag_co2_ice(:,:) = -1.
875   endif
876
877! Perennial CO2 ice
878    perennial_co2ice = 0.
879    call get_field("perennial_co2ice",perennial_co2ice,found,indextime)
880    if (.not. found) then
881        write(*,*) "phyetat0: Failed loading <perennial_co2ice> : ", &
882                   "<perennial_co2ice> is set as 10m at the South Pole"
883        if (abs(latitude(ngrid) - (-pi/2.)) < 1.e-5) perennial_co2ice(ngrid,:) = 10*1.6e3 ! 10m which is convert to kg/m^2
884    endif ! not found
885  else ! no startphy_file
886    h2o_ice_depth = -1.
887    lag_co2_ice = -1.
888    zdqsdif_ssi_tot = 0.
889    d_coef = 4.e-4
890    perennial_co2ice = 0.
891    if (abs(latitude(ngrid) - (-pi/2.)) < 1.e-5) perennial_co2ice(ngrid,:) = 10*1.6e3 ! 10m which is convert to kg/m^2
892  endif !startphy_file
893else
894   h2o_ice_depth = -1.
895   lag_co2_ice = -1.
896   zdqsdif_ssi_tot = 0.
897   d_coef = 4.e-4
898   perennial_co2ice = 0.
899endif !paleoclimate
900
901! close file:
902if (startphy_file) call close_startphy
903
904end subroutine phyetat0
905
906!======================================================================
907subroutine ini_tab_controle_dyn_xios(idayref)
908
909use comcstfi_h,        only: g, mugaz, omeg, rad, rcp
910use time_phylmdz_mod,  only: hour_ini, daysec, dtphys
911use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
912
913implicit none
914
915integer*4, intent(in) :: idayref ! date (initial date for this run)
916
917integer :: length, l
918parameter (length = 100)
919real    :: tab_cntrl(length) ! run parameters are stored in this array
920
921do l = 1,length
922    tab_cntrl(l) = 0.
923enddo
924tab_cntrl(1)  = real(nbp_lon)
925tab_cntrl(2)  = real(nbp_lat-1)
926tab_cntrl(3)  = real(nbp_lev)
927tab_cntrl(4)  = real(idayref)
928tab_cntrl(5)  = rad
929tab_cntrl(6)  = omeg
930tab_cntrl(7)  = g
931tab_cntrl(8)  = mugaz
932tab_cntrl(9)  = rcp
933tab_cntrl(10) = daysec
934tab_cntrl(11) = dtphys
935
936tab_cntrl(27) = hour_ini
937
938tab_cntrl_mod = tab_cntrl
939
940end subroutine ini_tab_controle_dyn_xios
941
942end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.