source: trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 @ 2613

Last change on this file since 2613 was 2557, checked in by emillour, 3 years ago

Generic GCM:
Fix in phyetat0, only call soil_settings if flag callsoil is true.
GM

File size: 13.3 KB
Line 
1module phyetat0_mod
2
3implicit none
4
5contains
6
7subroutine phyetat0 (startphy_file, &
8                     ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
9                     day_ini,time,tsurf,tsoil, &
10                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice, &
11                     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
12
13! to use  'getin_p'
14      use ioipsl_getin_p_mod, only: getin_p
15
16  use tabfi_mod, only: tabfi
17  USE tracer_h, ONLY: noms
18  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
19  use iostart, only: nid_start, open_startphy, close_startphy, &
20                     get_field, get_var, inquire_field, &
21                     inquire_dimension, inquire_dimension_length
22  use slab_ice_h, only: noceanmx
23  use callkeys_mod, only: CLFvarying,surfalbedo,surfemis, callsoil
24  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
25                                east_gwstress, west_gwstress
26  implicit none
27
28!======================================================================
29!  Arguments:
30!  ---------
31!  inputs:
32  logical,intent(in) :: startphy_file ! .true. if reading start file
33  integer,intent(in) :: ngrid
34  integer,intent(in) :: nlayer
35  character*(*),intent(in) :: fichnom ! "startfi.nc" file
36  integer,intent(in) :: tab0
37  integer,intent(in) :: Lmodif
38  integer,intent(in) :: nsoil ! # of soil layers
39  integer,intent(in) :: nq
40  integer,intent(out) :: day_ini
41  real,intent(out) :: time
42
43!  outputs:
44  real,intent(out) :: tsurf(ngrid) ! surface temperature
45  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
46  real,intent(out) :: emis(ngrid) ! surface emissivity
47  real,intent(out) :: q2(ngrid,nlayer+1) !
48  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
49! real co2ice(ngrid) ! co2 ice cover
50  real,intent(out) :: cloudfrac(ngrid,nlayer)
51  real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
52  real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,noceanmx) 
53  real,intent(out) :: tsea_ice(ngrid),sea_ice(ngrid)
54  real,intent(out) :: rnat(ngrid)
55
56!======================================================================
57!  Local variables:
58
59!      INTEGER radpas
60!      REAL co2_ppm
61!      REAL solaire
62
63      real xmin,xmax ! to display min and max of a field
64!
65      INTEGER ig,iq,lmax
66      INTEGER nid, nvarid
67      INTEGER ierr, i, nsrf
68!      integer isoil
69!      INTEGER length
70!      PARAMETER (length=100)
71      CHARACTER*7 str7
72      CHARACTER*2 str2
73      CHARACTER*1 yes
74!
75      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
76      INTEGER nqold
77
78! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
79!      logical :: oldtracernames=.false.
80      integer :: count
81      character(len=30) :: txt ! to store some text
82     
83      INTEGER :: indextime=1 ! index of selected time, default value=1
84      logical :: found
85     
86      character(len=8) :: modname="phyetat0"
87
88!
89! ALLOCATE ARRAYS IN surfdat_h
90!
91IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
92IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
93IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
94IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
95IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
96IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
97IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
98
99if (startphy_file) then
100  ! open physics initial state file:
101  call open_startphy(fichnom)
102
103  ! possibility to modify tab_cntrl in tabfi
104  write(*,*)
105  write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
106  call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
107                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
108
109else ! "academic" initialization of planetary parameters via tabfi
110  call tabfi (ngrid,0,0,0,day_ini,lmax,p_rad, &
111                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
112endif ! of if (startphy_file)
113
114if (startphy_file) then
115  ! Load surface geopotential:
116  call get_field("phisfi",phisfi,found)
117  if (.not.found) then
118    call abort_physic(modname,"Failed loading <phisfi>",1)
119  endif
120else
121  phisfi(:)=0.
122endif ! of if (startphy_file)
123write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
124               minval(phisfi), maxval(phisfi)
125
126if (startphy_file) then
127  ! Load bare ground albedo:
128  call get_field("albedodat",albedodat,found)
129  if (.not.found) then
130    call abort_physic(modname,"Failed loading <albedodat>",1)
131  endif
132else
133  ! If no startfi file, use parameter surfalbedo in def file
134  surfalbedo=0.5
135  call getin_p("surfalbedo",surfalbedo)
136  print*,"surfalbedo",surfalbedo
137  albedodat(:)=surfalbedo
138endif ! of if (startphy_file)
139write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
140             minval(albedodat), maxval(albedodat)
141
142! ZMEA
143if (startphy_file) then
144  call get_field("ZMEA",zmea,found)
145  if (.not.found) then
146    call abort_physic(modname,"Failed loading <ZMEA>",1)
147  endif
148else
149  zmea(:)=0.
150endif ! of if (startphy_file)
151write(*,*) "phyetat0: <ZMEA> range:", &
152             minval(zmea), maxval(zmea)
153
154! ZSTD
155if (startphy_file) then
156  call get_field("ZSTD",zstd,found)
157  if (.not.found) then
158    call abort_physic(modname,"Failed loading <ZSTD>",1)
159  endif
160else
161  zstd(:)=0.
162endif ! of if (startphy_file)
163write(*,*) "phyetat0: <ZSTD> range:", &
164             minval(zstd), maxval(zstd)
165
166! ZSIG
167if (startphy_file) then
168  call get_field("ZSIG",zsig,found)
169  if (.not.found) then
170    call abort_physic(modname,"Failed loading <ZSIG>",1)
171  endif
172else
173  zsig(:)=0.
174endif ! of if (startphy_file)
175write(*,*) "phyetat0: <ZSIG> range:", &
176             minval(zsig), maxval(zsig)
177
178! ZGAM
179if (startphy_file) then
180  call get_field("ZGAM",zgam,found)
181  if (.not.found) then
182    call abort_physic(modname,"Failed loading <ZGAM>",1)
183  endif
184else
185  zgam(:)=0.
186endif ! of if (startphy_file)
187write(*,*) "phyetat0: <ZGAM> range:", &
188             minval(zgam), maxval(zgam)
189
190! ZTHE
191if (startphy_file) then
192  call get_field("ZTHE",zthe,found)
193  if (.not.found) then
194    call abort_physic(modname,"Failed loading <ZTHE>",1)
195  endif
196else
197  zthe(:)=0.
198endif ! of if (startphy_file)
199write(*,*) "phyetat0: <ZTHE> range:", &
200             minval(zthe), maxval(zthe)
201
202! Surface temperature :
203if (startphy_file) then
204  call get_field("tsurf",tsurf,found,indextime)
205  if (.not.found) then
206    call abort_physic(modname,"Failed loading <tsurf>",1)
207  endif
208else
209  tsurf(:)=0. ! will be updated afterwards in physiq !
210endif ! of if (startphy_file)
211write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
212             minval(tsurf), maxval(tsurf)
213
214! Surface emissivity
215if (startphy_file) then
216  call get_field("emis",emis,found,indextime)
217  if (.not.found) then
218    call abort_physic(modname,"Failed loading <emis>",1)
219  endif
220else
221  ! If no startfi file, use parameter surfemis in def file
222  surfemis=1.0
223  call getin_p("surfemis",surfemis)
224  print*,"surfemis",surfemis
225  emis(:)=surfemis
226endif ! of if (startphy_file)
227write(*,*) "phyetat0: Surface emissivity <emis> range:", &
228             minval(emis), maxval(emis)
229
230! Cloud fraction (added by BC 2010)
231if (CLFvarying) then
232  if (startphy_file) then
233    call get_field("cloudfrac",cloudfrac,found,indextime)
234    if (.not.found) then
235      call abort_physic(modname,"Failed loading <cloudfrac>",1)
236    endif
237  else
238    cloudfrac(:,:)=0.0
239  endif ! of if (startphy_file)
240  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
241               minval(cloudfrac), maxval(cloudfrac)
242else
243  cloudfrac(:,:)=0.0
244endif ! of if (CLFvarying)
245
246! Total cloud fraction (added by BC 2010)
247if (CLFvarying) then
248  if (startphy_file) then
249    call get_field("totcloudfrac",totcloudfrac,found,indextime)
250    if (.not.found) then
251      call abort_physic(modname,"Failed loading <totcloudfrac>",1)
252    endif
253  else
254    totcloudfrac(:)=0.0
255  endif ! of if (startphy_file)
256  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
257             minval(totcloudfrac), maxval(totcloudfrac)
258else
259  totcloudfrac(:)=0.0
260endif ! of if (CLFvarying)
261
262! Height of oceanic ice (added by BC 2010)
263if (startphy_file) then
264  call get_field("hice",hice,found,indextime)
265  if (.not.found) then
266    write(*,*) "phyetat0: Failed loading <hice>"
267    !  call abort
268    hice(:)=0.
269  endif
270else
271  hice(:)=0.
272endif ! of if (startphy_file)
273write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
274             minval(hice), maxval(hice)
275
276! SLAB OCEAN (added by BC 2014)
277if (startphy_file) then
278  ! nature of the surface
279  call get_field("rnat",rnat,found,indextime)
280  if (.not.found) then
281    write(*,*) "phyetat0: Failed loading <rnat>"
282        rnat(1:ngrid)=1.
283  else
284      do ig=1,ngrid
285        if((nint(rnat(ig)).eq.2).or.(nint(rnat(ig)).eq.0))then
286          rnat(ig)=0.
287        else
288          rnat(ig)=1.
289        endif     
290      enddo
291  endif ! of if (.not.found)
292else
293  rnat(:)=1.
294endif ! of if (startphy_file)
295write(*,*) "phyetat0: Nature of surface <rnat> range:", &
296             minval(rnat), maxval(rnat)
297
298if (startphy_file) then
299  ! Pourcentage of sea ice cover
300  call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
301  if (.not.found) then
302    write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
303    pctsrf_sic(1:ngrid)=0.
304  endif
305else
306  pctsrf_sic(:)=0.
307endif ! of if (startphy_file)
308write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
309             minval(pctsrf_sic), maxval(pctsrf_sic)
310
311if (startphy_file) then
312  ! Slab ocean temperature (2 layers)
313  call get_field("tslab",tslab,found,indextime)
314  if (.not.found) then
315    write(*,*) "phyetat0: Failed loading <tslab>"
316    do iq=1,noceanmx
317      tslab(1:ngrid,iq)=tsurf(1:ngrid)
318    enddo
319  endif
320else
321  do iq=1,noceanmx
322    tslab(1:ngrid,iq)=tsurf(1:ngrid)
323  enddo
324endif ! of if (startphy_file)
325write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
326             minval(tslab), maxval(tslab)
327
328if (startphy_file) then
329  ! Oceanic ice temperature
330  call get_field("tsea_ice",tsea_ice,found,indextime)
331  if (.not.found) then
332    write(*,*) "phyetat0: Failed loading <tsea_ice>"
333    tsea_ice(1:ngrid)=273.15-1.8
334  endif
335else
336  tsea_ice(1:ngrid)=273.15-1.8
337endif ! of if (startphy_file)
338write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
339             minval(tsea_ice), maxval(tsea_ice)
340
341if (startphy_file) then
342  !  Oceanic ice quantity (kg/m^2)
343  call get_field("sea_ice",sea_ice,found,indextime)
344  if (.not.found) then
345    write(*,*) "phyetat0: Failed loading <sea_ice>"
346    sea_ice(1:ngrid)=0.
347  endif
348else
349  sea_ice(1:ngrid)=0.
350endif ! of if (startphy_file)
351write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
352             minval(sea_ice), maxval(sea_ice)
353
354
355! pbl wind variance
356if (startphy_file) then
357  call get_field("q2",q2,found,indextime)
358  if (.not.found) then
359    call abort_physic(modname,"Failed loading <q2>",1)
360  endif
361else
362  q2(:,:)=0.
363endif ! of if (startphy_file)
364write(*,*) "phyetat0: PBL wind variance <q2> range:", &
365             minval(q2), maxval(q2)
366
367! tracer on surface
368if (nq.ge.1) then
369  do iq=1,nq
370    txt=noms(iq)
371    if (startphy_file) then
372      call get_field(txt,qsurf(:,iq),found,indextime)
373      if (.not.found) then
374        write(*,*) "phyetat0: Failed loading <",trim(txt),">"
375        write(*,*) "         ",trim(txt)," is set to zero"
376        qsurf(:,iq) = 0.
377      endif
378    else
379      qsurf(:,iq)=0.
380    endif ! of if (startphy_file)
381    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
382                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
383  enddo! of do iq=1,nq
384endif ! of if (nq.ge.1)
385
386
387if (startphy_file) then
388  ! Call to soil_settings, in order to read soil temperatures,
389  ! as well as thermal inertia and volumetric heat capacity
390  if (callsoil) then
391    call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
392  endif
393endif ! of if (startphy_file)
394
395! Non-orographic gravity waves
396if (startphy_file) then
397   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
398   if (.not.found) then
399      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
400      du_nonoro_gwd(:,:)=0.
401   else
402      write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
403      write(*,*) " <du_nonoro_gwd> range:", &
404             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
405   endif
406endif ! of if (startphy_file)
407if (startphy_file) then
408   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
409   if (.not.found) then
410      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
411      dv_nonoro_gwd(:,:)=0.
412   else
413      write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
414      write(*,*) " <dv_nonoro_gwd> range:", &
415             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
416   endif
417endif ! of if (startphy_file)
418if (startphy_file) then
419   call get_field("east_gwstress",east_gwstress,found,indextime)
420   if (.not.found) then
421      write(*,*) "phyetat0: <east_gwstress> not in file"
422      east_gwstress(:,:)=0.
423   else
424      write(*,*) "phyetat0: Memory of eastward stress profile due to non-orographic GW"
425      write(*,*) " <east_gwstress> range:", &
426             minval(east_gwstress), maxval(east_gwstress)
427   endif
428endif ! of if (startphy_file)
429if (startphy_file) then
430   call get_field("west_gwstress",west_gwstress,found,indextime)
431   if (.not.found) then
432      write(*,*) "phyetat0: <west_gwstress> not in file"
433      west_gwstress(:,:)=0.
434   else
435      write(*,*) "phyetat0: Memory of westward stress profile due to non-orographic GW"
436      write(*,*) " <west_gwstress> range:", &
437             minval(west_gwstress), maxval(west_gwstress)
438   endif
439endif ! of if (startphy_file)
440!
441! close file:
442!
443if (startphy_file) call close_startphy
444
445end subroutine phyetat0
446
447end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.