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

Last change on this file since 2220 was 2220, checked in by dbardet, 5 years ago

Update the nonorographic gravity waves drag parametrization (from the r3599 of Earth s model LMDZ6): 1) add east_ and west_gwdstress variables, 2) delete aleas function: random waves are producted by the MOD function, 3) reproductibility concerning the number of procs is now validated, 4) changing name of some variables like RUW, RVW, ZOP, ZOM,..., 5) tendency of winds due to GW drag are now module variables and written in restart files.
File size: 13.8 KB
Line 
1module phyetat0_mod
2
3implicit none
4
5contains
6
7subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
8                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf,co2ice, &
9                     tauscaling,totcloudfrac,wstar,mem_Mccn_co2,mem_Nccn_co2,&
10                     mem_Mh2o_co2)
11
12  use tracer_mod, only: noms ! tracer names
13  use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
14                       zmea, zstd, zsig, zgam, zthe, hmons, summit, base
15  use iostart, only: nid_start, open_startphy, close_startphy, &
16                     get_field, get_var, inquire_field, &
17                     inquire_dimension, inquire_dimension_length
18  use ioipsl_getincom, only : getin
19  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
20  implicit none
21 
22  include "callkeys.h"
23!======================================================================
24! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
25!  Adaptation à Mars : Yann Wanherdrick
26! Objet: Lecture de l etat initial pour la physique
27! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
28!                      r4 or r8 restarts independently of having compiled
29!                      the GCM in r4 or r8)
30!         June 2013 TN : Possibility to read files with a time axis
31!         November 2013 EM : Enabeling parallel, using iostart module
32!======================================================================
33  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
34  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
35!======================================================================
36!  Arguments:
37!  ---------
38!  inputs:
39  character*(*),intent(in) :: fichnom ! "startfi.nc" file
40  integer,intent(in) :: tab0
41  integer,intent(in) :: Lmodif
42  integer,intent(in) :: nsoil ! # of soil layers
43  integer,intent(in) :: ngrid ! # of atmospheric columns
44  integer,intent(in) :: nlay ! # of atmospheric layers
45  integer,intent(in) :: nq
46  integer :: day_ini
47  real :: time0
48
49!  outputs:
50  real,intent(out) :: tsurf(ngrid) ! surface temperature
51  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
52  real,intent(out) :: albedo(ngrid,2) ! surface albedo
53  real,intent(out) :: emis(ngrid) ! surface emissivity
54  real,intent(out) :: q2(ngrid,nlay+1) !
55  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
56  real,intent(out) :: co2ice(ngrid) ! co2 ice cover
57  real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
58  real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
59  real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
60  real,intent(out) :: mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
61  real,intent(out) :: mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
62  real,intent(out) :: mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
63!======================================================================
64!  Local variables:
65
66      real surffield(ngrid) ! to temporarily store a surface field
67      real xmin,xmax ! to display min and max of a field
68!
69      INTEGER ig,iq,lmax
70      INTEGER nid, nvarid
71      INTEGER ierr, i, nsrf
72!      integer isoil
73!      INTEGER length
74!      PARAMETER (length=100)
75      CHARACTER*7 str7
76      CHARACTER*2 str2
77      CHARACTER*1 yes
78!
79      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
80      INTEGER nqold
81
82! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
83      logical :: oldtracernames=.false.
84      integer :: count
85      character(len=30) :: txt ! to store some text
86
87! specific for time
88      REAL,ALLOCATABLE :: time(:) ! times stored in start
89      INTEGER timelen ! number of times stored in the file
90      INTEGER indextime ! index of selected time
91
92      INTEGER :: edges(3),corner(3)
93      LOGICAL :: found
94
95      REAL :: timestart ! to pick which initial state to start from
96
97! open physics initial state file:
98call open_startphy(fichnom)
99
100
101! possibility to modify tab_cntrl in tabfi
102write(*,*)
103write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
104call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
105            p_omeg,p_g,p_mugaz,p_daysec,time0)
106
107
108! Load surface geopotential:
109call get_field("phisfi",phisfi,found)
110if (.not.found) then
111  write(*,*) "phyetat0: Failed loading <phisfi>"
112  call abort
113else
114  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
115             minval(phisfi), maxval(phisfi)
116endif
117
118
119! Load bare ground albedo:
120call get_field("albedodat",albedodat,found)
121if (.not.found) then
122  write(*,*) "phyetat0: Failed loading <albedodat>"
123  call abort
124else
125  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
126             minval(albedodat), maxval(albedodat)
127endif
128
129! ZMEA
130call get_field("ZMEA",zmea,found)
131if (.not.found) then
132  write(*,*) "phyetat0: Failed loading <ZMEA>"
133  call abort
134else
135  write(*,*) "phyetat0: <ZMEA> range:", &
136             minval(zmea), maxval(zmea)
137endif
138
139
140! ZSTD
141call get_field("ZSTD",zstd,found)
142if (.not.found) then
143  write(*,*) "phyetat0: Failed loading <ZSTD>"
144  call abort
145else
146  write(*,*) "phyetat0: <ZSTD> range:", &
147             minval(zstd), maxval(zstd)
148endif
149
150
151! ZSIG
152call get_field("ZSIG",zsig,found)
153if (.not.found) then
154  write(*,*) "phyetat0: Failed loading <ZSIG>"
155  call abort
156else
157  write(*,*) "phyetat0: <ZSIG> range:", &
158             minval(zsig), maxval(zsig)
159endif
160
161
162! ZGAM
163call get_field("ZGAM",zgam,found)
164if (.not.found) then
165  write(*,*) "phyetat0: Failed loading <ZGAM>"
166  call abort
167else
168  write(*,*) "phyetat0: <ZGAM> range:", &
169             minval(zgam), maxval(zgam)
170endif
171
172
173! ZTHE
174call get_field("ZTHE",zthe,found)
175if (.not.found) then
176  write(*,*) "phyetat0: Failed loading <ZTHE>"
177  call abort
178else
179  write(*,*) "phyetat0: <ZTHE> range:", &
180             minval(zthe), maxval(zthe)
181endif
182
183! hmons
184call get_field("hmons",hmons,found)
185if (.not.found) then
186  write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
187  if (rdstorm) then
188    call abort
189  else
190    write(*,*) "will continue anyway..."
191    write(*,*) "because you may not need it."
192    hmons(:)=0.
193  end if
194else
195  do ig=1,ngrid
196    if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
197  enddo
198  write(*,*) "phyetat0: <hmons> range:", &
199             minval(hmons), maxval(hmons)
200endif
201
202! summit
203call get_field("summit",summit,found)
204if (.not.found) then
205  write(*,*) "WARNING: phyetat0: Failed loading <summit>"
206  if (rdstorm) then
207    call abort
208  else
209    write(*,*) "will continue anyway..."
210    write(*,*) "because you may not need it."
211    summit(:)=0.
212  end if
213else
214  do ig=1,ngrid
215    if (summit(ig) .eq. -999999.)  summit(ig)=0.
216  enddo
217  write(*,*) "phyetat0: <summit> range:", &
218             minval(summit), maxval(summit)
219endif
220
221! summit
222call get_field("base",base,found)
223if (.not.found) then
224  write(*,*) "WARNING: phyetat0: Failed loading <base>"
225  if (rdstorm) then
226    call abort
227  else
228    write(*,*) "will continue anyway..."
229    write(*,*) "because you may not need it."
230    base(:)=0.
231  end if
232else
233  do ig=1,ngrid
234    if (base(ig) .eq. -999999.)  base(ig)=0.
235  enddo
236  write(*,*) "phyetat0: <base> range:", &
237             minval(base), maxval(base)
238endif
239 
240! Time axis
241! obtain timestart from run.def
242timestart=-9999 ! default value
243call getin("timestart",timestart)
244
245found=inquire_dimension("Time")
246if (.not.found) then
247  indextime = 1
248  write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
249else
250  write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
251  timelen=inquire_dimension_length("Time")
252  allocate(time(timelen))
253  ! load "Time" array:
254  call get_var("Time",time,found)
255  if (.not.found) then
256    write(*,*) "phyetat0: Failed loading <Time>"
257    call abort
258  endif
259  ! seclect the desired time index
260  IF (timestart .lt. 0) THEN  ! default: we use the last time value
261    indextime = timelen
262  ELSE  ! else we look for the desired value in the time axis
263    indextime = 0
264    DO i=1,timelen
265      IF (abs(time(i) - timestart) .lt. 0.01) THEN
266        indextime = i
267        EXIT
268      ENDIF
269    ENDDO
270    IF (indextime .eq. 0) THEN
271      PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
272      PRINT*, "Stored times are:"
273      DO i=1,timelen
274         PRINT*, time(i)
275      ENDDO
276      CALL abort
277    ENDIF
278  ENDIF ! of IF (timestart .lt. 0)
279  ! In startfi the absolute date is day_ini + time0 + time
280  ! For now on, in the GCM physics, it is day_ini + time0
281  time0 = time(indextime) + time0
282  day_ini = day_ini + INT(time0)
283  time0 = time0 - INT(time0)
284       
285  PRINT*, "phyetat0: Selected time ",time(indextime), &
286          " at index ",indextime
287     
288  DEALLOCATE(time)
289endif ! of if Time not found in file
290
291
292! CO2 ice cover
293call get_field("co2ice",co2ice,found,indextime)
294if (.not.found) then
295  write(*,*) "phyetat0: Failed loading <co2ice>"
296  call abort
297else
298  write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
299             minval(co2ice), maxval(co2ice)
300endif
301
302! Memory of the origin of the co2 particles
303call get_field("mem_Mccn_co2",mem_Mccn_co2,found,indextime)
304if (.not.found) then
305  write(*,*) "phyetat0: <mem_Mccn_co2> not in file"
306  mem_Mccn_co2(:,:)=0
307else
308  write(*,*) "phyetat0: Memory of CCN mass of H2O and dust used by CO2"
309  write(*,*) " <mem_Mccn_co2> range:", &
310             minval(mem_Mccn_co2), maxval(mem_Mccn_co2)
311endif
312
313call get_field("mem_Nccn_co2",mem_Nccn_co2,found,indextime)
314if (.not.found) then
315  write(*,*) "phyetat0: <mem_Nccn_co2> not in file"
316  mem_Nccn_co2(:,:)=0
317else
318  write(*,*) "phyetat0: Memory of CCN number of H2O and dust used by CO2"
319  write(*,*) " <mem_Nccn_co2> range:", &
320             minval(mem_Nccn_co2), maxval(mem_Nccn_co2)
321endif
322
323call get_field("mem_Mh2o_co2",mem_Mh2o_co2,found,indextime)
324if (.not.found) then
325  write(*,*) "phyetat0: <mem_Mh2o_co2> not in file"
326  mem_Mh2o_co2(:,:)=0
327else
328  write(*,*) "phyetat0: Memory of H2O mass integred into CO2 crystal"
329  write(*,*) " <mem_Mh2o_co2> range:", &
330             minval(mem_Mh2o_co2), maxval(mem_Mh2o_co2)
331endif
332
333! Dust conversion factor
334call get_field("tauscaling",tauscaling,found,indextime)
335if (.not.found) then
336  write(*,*) "phyetat0: <tauscaling> not in file"
337  tauscaling(:) = 1
338else
339  write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
340             minval(tauscaling), maxval(tauscaling)
341endif
342
343! Sub-grid cloud fraction
344call get_field("totcloudfrac",totcloudfrac,found,indextime)
345if (.not.found) then
346  write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
347  totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
348else
349  write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
350             minval(totcloudfrac), maxval(totcloudfrac)
351endif
352
353! Max vertical velocity in thermals
354call get_field("wstar",wstar,found,indextime)
355if (.not.found) then
356  write(*,*) "phyetat0: <wstar> not in file! Set to zero"
357  wstar(:)=0
358else
359  write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
360             minval(wstar),maxval(wstar)
361endif
362
363! Surface temperature :
364call get_field("tsurf",tsurf,found,indextime)
365if (.not.found) then
366  write(*,*) "phyetat0: Failed loading <tsurf>"
367  call abort
368else
369  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
370             minval(tsurf), maxval(tsurf)
371endif
372
373! Surface albedo
374call get_field("albedo",albedo(:,1),found,indextime)
375if (.not.found) then
376  write(*,*) "phyetat0: Failed loading <albedo>"
377  albedo(:,1)=albedodat(:)
378else
379  write(*,*) "phyetat0: Surface albedo <albedo> range:", &
380             minval(albedo(:,1)), maxval(albedo(:,1))
381endif
382albedo(:,2)=albedo(:,1)
383
384! Surface emissivity
385call get_field("emis",emis,found,indextime)
386if (.not.found) then
387  write(*,*) "phyetat0: Failed loading <emis>"
388  call abort
389else
390  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
391             minval(emis), maxval(emis)
392endif
393
394
395! surface roughness length (NB: z0 is a common in surfdat_h)
396call get_field("z0",z0,found)
397if (.not.found) then
398  write(*,*) "phyetat0: Failed loading <z0>"
399  write(*,*) 'will use constant value of z0_default:',z0_default
400  z0(:)=z0_default
401else
402  write(*,*) "phyetat0: Surface roughness <z0> range:", &
403             minval(z0), maxval(z0)
404endif
405
406
407! pbl wind variance
408call get_field("q2",q2,found,indextime)
409if (.not.found) then
410  write(*,*) "phyetat0: Failed loading <q2>"
411  call abort
412else
413  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
414             minval(q2), maxval(q2)
415endif
416
417! Non-orographic gravity waves
418call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
419if (.not.found) then
420   write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
421   du_nonoro_gwd(:,:)=0.
422else
423   write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
424   write(*,*) " <du_nonoro_gwd> range:", &
425          minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
426endif
427
428call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
429if (.not.found) then
430   write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
431   dv_nonoro_gwd(:,:)=0.
432else
433   write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
434   write(*,*) " <dv_nonoro_gwd> range:", &
435          minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
436endif
437
438! tracer on surface
439if (nq.ge.1) then
440  do iq=1,nq
441    txt=noms(iq)
442    if (txt.eq."h2o_vap") then
443      ! There is no surface tracer for h2o_vap;
444      ! "h2o_ice" should be loaded instead
445      txt="h2o_ice"
446      write(*,*) 'phyetat0: loading surface tracer', &
447                           ' h2o_ice instead of h2o_vap'
448    endif
449    call get_field(txt,qsurf(:,iq),found,indextime)
450    if (.not.found) then
451      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
452      write(*,*) "         ",trim(txt)," is set to zero"
453    else
454      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
455                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
456    endif
457  enddo
458endif ! of if (nq.ge.1)
459
460! Call to soil_settings, in order to read soil temperatures,
461! as well as thermal inertia and volumetric heat capacity
462call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
463
464!
465! close file:
466!
467call close_startphy
468
469end subroutine phyetat0
470
471end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.