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

Last change on this file since 3552 was 3552, checked in by emillour, 3 weeks ago

Generic PCM

Modify iostart routines to take netcdf index as an argument in order to
be used to read/write in other files than startfi.nc (preparing 1D
restart)

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