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

Last change on this file since 3533 was 3522, checked in by afalco, 2 months ago

Generic PCM: imported write_output() subroutine from Mars PCM to write both diagfi and XIOS.

Still some bug on a few variables from generic condensation, to be checked.

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