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

Last change on this file since 3393 was 3393, checked in by jliu, 5 months ago

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