source: trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F90 @ 1477

Last change on this file since 1477 was 1421, checked in by milmd, 10 years ago

Correction of newstart and vertical interpolation of q2 in LMDZ.GENERIC and LMDZ.MARS. In LMDZ.GENERIC add flag for cloud fraction reading.
File size: 11.4 KB
Line 
1subroutine phyetat0 (ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
2                     day_ini,time,tsurf,tsoil, &
3                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice, &
4                     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
5
6
7  USE infotrac, ONLY: tname
8  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
9  use iostart, only: nid_start, open_startphy, close_startphy, &
10                     get_field, get_var, inquire_field, &
11                     inquire_dimension, inquire_dimension_length
12  use slab_ice_h, only: noceanmx
13  use callkeys_mod, only: CLFvarying
14
15  implicit none
16
17!======================================================================
18! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
19!  Adaptation à Mars : Yann Wanherdrick
20! Objet: Lecture de l etat initial pour la physique
21!======================================================================
22#include "netcdf.inc"
23
24!======================================================================
25!  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
26!  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
27!======================================================================
28!  Arguments:
29!  ---------
30!  inputs:
31  integer,intent(in) :: ngrid
32  integer,intent(in) :: nlayer
33  character*(*),intent(in) :: fichnom ! "startfi.nc" file
34  integer,intent(in) :: tab0
35  integer,intent(in) :: Lmodif
36  integer,intent(in) :: nsoil ! # of soil layers
37  integer,intent(in) :: nq
38  integer,intent(in) :: day_ini
39  real,intent(in) :: time
40
41!  outputs:
42  real,intent(out) :: tsurf(ngrid) ! surface temperature
43  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
44  real,intent(out) :: emis(ngrid) ! surface emissivity
45  real,intent(out) :: q2(ngrid,nlayer+1) !
46  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
47! real co2ice(ngrid) ! co2 ice cover
48  real,intent(out) :: cloudfrac(ngrid,nlayer)
49  real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
50  real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,noceanmx) 
51  real,intent(out) :: tsea_ice(ngrid),sea_ice(ngrid)
52  real,intent(out) :: rnat(ngrid)
53
54!======================================================================
55!  Local variables:
56
57!      INTEGER radpas
58!      REAL co2_ppm
59!      REAL solaire
60
61      real xmin,xmax ! to display min and max of a field
62!
63      INTEGER ig,iq,lmax
64      INTEGER nid, nvarid
65      INTEGER ierr, i, nsrf
66!      integer isoil
67!      INTEGER length
68!      PARAMETER (length=100)
69      CHARACTER*7 str7
70      CHARACTER*2 str2
71      CHARACTER*1 yes
72!
73      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
74      INTEGER nqold
75
76! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
77!      logical :: oldtracernames=.false.
78      integer :: count
79      character(len=30) :: txt ! to store some text
80     
81      INTEGER :: indextime=1 ! index of selected time, default value=1
82      logical :: found
83
84!
85! ALLOCATE ARRAYS IN surfdat_h
86!
87IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
88IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
89IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
90IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
91IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
92IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
93IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
94
95
96! open physics initial state file:
97call open_startphy(fichnom)
98
99
100! possibility to modify tab_cntrl in tabfi
101write(*,*)
102write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
103call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
104                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
105
106!c
107!c Lecture des latitudes (coordonnees):
108!c
109!      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
110!      IF (ierr.NE.NF_NOERR) THEN
111!         PRINT*, 'phyetat0: Le champ <latitude> est absent'
112!         CALL abort
113!      ENDIF
114!#ifdef NC_DOUBLE
115!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
116!#else
117!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
118!#endif
119!      IF (ierr.NE.NF_NOERR) THEN
120!         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
121!         CALL abort
122!      ENDIF
123!c
124!c Lecture des longitudes (coordonnees):
125!c
126!      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
127!      IF (ierr.NE.NF_NOERR) THEN
128!         PRINT*, 'phyetat0: Le champ <longitude> est absent'
129!         CALL abort
130!      ENDIF
131!#ifdef NC_DOUBLE
132!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
133!#else
134!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
135!#endif
136!      IF (ierr.NE.NF_NOERR) THEN
137!         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
138!         CALL abort
139!      ENDIF
140!c
141!c Lecture des aires des mailles:
142!c
143!      ierr = NF_INQ_VARID (nid, "area", nvarid)
144!      IF (ierr.NE.NF_NOERR) THEN
145!         PRINT*, 'phyetat0: Le champ <area> est absent'
146!         CALL abort
147!      ENDIF
148!#ifdef NC_DOUBLE
149!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
150!#else
151!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
152!#endif
153!      IF (ierr.NE.NF_NOERR) THEN
154!         PRINT*, 'phyetat0: Lecture echouee pour <area>'
155!         CALL abort
156!      ENDIF
157!      xmin = 1.0E+20
158!      xmax = -1.0E+20
159!      xmin = MINVAL(area)
160!      xmax = MAXVAL(area)
161!      PRINT*,'Aires des mailles <area>:', xmin, xmax
162
163! Load surface geopotential:
164call get_field("phisfi",phisfi,found)
165if (.not.found) then
166  write(*,*) "phyetat0: Failed loading <phisfi>"
167  call abort
168else
169  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
170             minval(phisfi), maxval(phisfi)
171endif
172
173! Load bare ground albedo:
174call get_field("albedodat",albedodat,found)
175if (.not.found) then
176  write(*,*) "phyetat0: Failed loading <albedodat>"
177  call abort
178else
179  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
180             minval(albedodat), maxval(albedodat)
181endif
182
183! ZMEA
184call get_field("ZMEA",zmea,found)
185if (.not.found) then
186  write(*,*) "phyetat0: Failed loading <ZMEA>"
187  call abort
188else
189  write(*,*) "phyetat0: <ZMEA> range:", &
190             minval(zmea), maxval(zmea)
191endif
192
193! ZSTD
194call get_field("ZSTD",zstd,found)
195if (.not.found) then
196  write(*,*) "phyetat0: Failed loading <ZSTD>"
197  call abort
198else
199  write(*,*) "phyetat0: <ZSTD> range:", &
200             minval(zstd), maxval(zstd)
201endif
202
203! ZSIG
204call get_field("ZSIG",zsig,found)
205if (.not.found) then
206  write(*,*) "phyetat0: Failed loading <ZSIG>"
207  call abort
208else
209  write(*,*) "phyetat0: <ZSIG> range:", &
210             minval(zsig), maxval(zsig)
211endif
212
213! ZGAM
214call get_field("ZGAM",zgam,found)
215if (.not.found) then
216  write(*,*) "phyetat0: Failed loading <ZGAM>"
217  call abort
218else
219  write(*,*) "phyetat0: <ZGAM> range:", &
220             minval(zgam), maxval(zgam)
221endif
222
223! ZTHE
224call get_field("ZTHE",zthe,found)
225if (.not.found) then
226  write(*,*) "phyetat0: Failed loading <ZTHE>"
227  call abort
228else
229  write(*,*) "phyetat0: <ZTHE> range:", &
230             minval(zthe), maxval(zthe)
231endif
232
233! Surface temperature :
234call get_field("tsurf",tsurf,found,indextime)
235if (.not.found) then
236  write(*,*) "phyetat0: Failed loading <tsurf>"
237  call abort
238else
239  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
240             minval(tsurf), maxval(tsurf)
241endif
242
243! Surface emissivity
244call get_field("emis",emis,found,indextime)
245if (.not.found) then
246  write(*,*) "phyetat0: Failed loading <emis>"
247  call abort
248else
249  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
250             minval(emis), maxval(emis)
251endif
252
253! Cloud fraction (added by BC 2010)
254if (CLFvarying) then
255call get_field("cloudfrac",cloudfrac,found,indextime)
256if (.not.found) then
257  write(*,*) "phyetat0: Failed loading <cloudfrac>"
258  call abort
259else
260  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
261             minval(cloudfrac), maxval(cloudfrac)
262endif
263else
264cloudfrac(:,:)=0.0
265endif
266
267! Total cloud fraction (added by BC 2010)
268if (CLFvarying) then
269call get_field("totcloudfrac",totcloudfrac,found,indextime)
270if (.not.found) then
271  write(*,*) "phyetat0: Failed loading <totcloudfrac>"
272  call abort
273else
274  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
275             minval(totcloudfrac), maxval(totcloudfrac)
276endif
277else
278totcloudfrac(:)=0.0
279endif
280
281! Height of oceanic ice (added by BC 2010)
282call get_field("hice",hice,found,indextime)
283if (.not.found) then
284  write(*,*) "phyetat0: Failed loading <hice>"
285!  call abort
286      do ig=1,ngrid
287      hice(ig)=0.
288      enddo
289else
290  write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
291             minval(hice), maxval(hice)
292endif
293
294! SLAB OCEAN (added by BC 2014)
295! nature of the surface
296call get_field("rnat",rnat,found,indextime)
297if (.not.found) then
298  write(*,*) "phyetat0: Failed loading <rnat>"
299      do ig=1,ngrid
300        rnat(ig)=1.
301      enddo
302else
303      do ig=1,ngrid
304        if((nint(rnat(ig)).eq.2).or.(nint(rnat(ig)).eq.0))then
305          rnat(ig)=0.
306        else
307          rnat(ig)=1.
308        endif     
309      enddo
310
311  write(*,*) "phyetat0: Nature of surface <rnat> range:", &
312             minval(rnat), maxval(rnat)
313endif
314! Pourcentage of sea ice cover
315call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
316if (.not.found) then
317  write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
318      do ig=1,ngrid
319      pctsrf_sic(ig)=0.
320      enddo
321else
322  write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
323             minval(pctsrf_sic), maxval(pctsrf_sic)
324endif
325! Slab ocean temperature (2 layers)
326call get_field("tslab",tslab,found,indextime)
327if (.not.found) then
328  write(*,*) "phyetat0: Failed loading <tslab>"
329      do ig=1,ngrid
330      do iq=1,noceanmx
331      tslab(ig,iq)=tsurf(ig)
332      enddo
333      enddo
334else
335  write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
336             minval(tslab), maxval(tslab)
337endif
338! Oceanic ice temperature
339call get_field("tsea_ice",tsea_ice,found,indextime)
340if (.not.found) then
341  write(*,*) "phyetat0: Failed loading <tsea_ice>"
342      do ig=1,ngrid
343      tsea_ice(ig)=273.15-1.8
344      enddo
345else
346  write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
347             minval(tsea_ice), maxval(tsea_ice)
348endif
349!  Oceanic ice quantity (kg/m^2)
350call get_field("sea_ice",sea_ice,found,indextime)
351if (.not.found) then
352  write(*,*) "phyetat0: Failed loading <sea_ice>"
353      do ig=1,ngrid
354      tsea_ice(ig)=0.
355      enddo
356else
357  write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
358             minval(sea_ice), maxval(sea_ice)
359endif
360
361
362
363
364! pbl wind variance
365call get_field("q2",q2,found,indextime)
366if (.not.found) then
367  write(*,*) "phyetat0: Failed loading <q2>"
368  call abort
369else
370  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
371             minval(q2), maxval(q2)
372endif
373
374! tracer on surface
375if (nq.ge.1) then
376  do iq=1,nq
377    txt=tname(iq)
378    if (txt.eq."h2o_vap") then
379      ! There is no surface tracer for h2o_vap;
380      ! "h2o_ice" should be loaded instead
381      txt="h2o_ice"
382      write(*,*) 'phyetat0: loading surface tracer', &
383                           ' h2o_ice instead of h2o_vap'
384    endif
385    call get_field(txt,qsurf(:,iq),found,indextime)
386    if (.not.found) then
387      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
388      write(*,*) "         ",trim(txt)," is set to zero"
389      qsurf(:,iq) = 0.
390    else
391      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
392                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
393    endif
394  enddo
395endif ! of if (nq.ge.1)
396
397
398! Call to soil_settings, in order to read soil temperatures,
399! as well as thermal inertia and volumetric heat capacity
400call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
401!
402! close file:
403!
404call close_startphy
405
406END SUBROUTINE phyetat0
Note: See TracBrowser for help on using the repository browser.