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

Last change on this file since 1957 was 1763, checked in by jvatant, 8 years ago

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