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

Last change on this file since 1962 was 1944, checked in by emillour, 7 years ago

Mars GCM:
A first step towards 1+1=2 (for now only works without tracers):

  • store and load "albedo" (surface albedo) and wstar (thermals' max vertical velocity) in physics (re)start file.
  • turn phyetat0 into a module in the process.

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