source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0.F90

Last change on this file was 313, checked in by ymipsl, 10 years ago
  • implement splitting of XIOS file for lmdz physics
  • Termination is done properly in parallel by calling MPI_ABORT instead of abort or stop

YM

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