source: trunk/LMDZ.MARS/libf/phymars/phyetat0.F90 @ 1930

Last change on this file since 1930 was 1922, checked in by emillour, 8 years ago

Mars GCM:
CO2 code updates:

  • make co2cloud a module and save mem_* variables (initialized via phys_state_var_init)
  • make improvedCO2cloud a module
  • read/write mem_* variables in phyetat0.F and phyredem.F

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