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

Last change on this file since 2266 was 2265, checked in by emillour, 5 years ago

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