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

Last change on this file since 3347 was 3335, checked in by emillour, 23 months ago

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