source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0_academic.F90 @ 232

Last change on this file since 232 was 232, checked in by ymipsl, 10 years ago

Work around for init etat0 state => missing field in the startfi are initialize to constant value.

YM

File size: 11.4 KB
Line 
1subroutine phyetat0_academic (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
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
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
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
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
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
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  phisfi(:)=0
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  do ig=1,ngrid
181    albedodat(ig)=0.
182  enddo
183else
184  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
185             minval(albedodat), maxval(albedodat)
186endif
187
188! ZMEA
189call get_field("ZMEA",zmea,found)
190if (.not.found) then
191    zmea(:)=0.
192else
193  write(*,*) "phyetat0: <ZMEA> range:", &
194             minval(zmea), maxval(zmea)
195endif
196
197! ZSTD
198call get_field("ZSTD",zstd,found)
199if (.not.found) then
200  write(*,*) "phyetat0: Failed loading <ZSTD>"
201   zstd(:)=0.
202else
203  write(*,*) "phyetat0: <ZSTD> range:", &
204             minval(zstd), maxval(zstd)
205endif
206
207! ZSIG
208call get_field("ZSIG",zsig,found)
209if (.not.found) then
210  write(*,*) "phyetat0: Failed loading <ZSIG>"
211  zsig(:)=0.
212else
213  write(*,*) "phyetat0: <ZSIG> range:", &
214             minval(zsig), maxval(zsig)
215endif
216
217! ZGAM
218call get_field("ZGAM",zgam,found)
219if (.not.found) then
220  write(*,*) "phyetat0: Failed loading <ZGAM>"
221  zgam(:)=0.
222else
223  write(*,*) "phyetat0: <ZGAM> range:", &
224             minval(zgam), maxval(zgam)
225endif
226
227! ZTHE
228call get_field("ZTHE",zthe,found)
229if (.not.found) then
230  write(*,*) "phyetat0: Failed loading <ZTHE>"
231  zthe(:)=0.
232else
233  write(*,*) "phyetat0: <ZTHE> range:", &
234             minval(zthe), maxval(zthe)
235endif
236
237! Surface temperature :
238call get_field("tsurf",tsurf,found,indextime)
239if (.not.found) then
240  write(*,*) "phyetat0: Failed loading <tsurf>"
241  tsurf(:)=175.0
242else
243  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
244             minval(tsurf), maxval(tsurf)
245endif
246
247! Surface emissivity
248call get_field("emis",emis,found,indextime)
249if (.not.found) then
250  write(*,*) "phyetat0: Failed loading <emis>"
251  emis(:)=0.5
252else
253  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
254             minval(emis), maxval(emis)
255endif
256
257! Cloud fraction (added by BC 2010)
258call get_field("cloudfrac",cloudfrac,found,indextime)
259if (.not.found) then
260  write(*,*) "phyetat0: Failed loading <cloudfrac>"
261  cloudfrac(:,:)=0.
262else
263  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
264             minval(cloudfrac), maxval(cloudfrac)
265endif
266
267! Total cloud fraction (added by BC 2010)
268call get_field("totcloudfrac",totcloudfrac,found,indextime)
269if (.not.found) then
270  write(*,*) "phyetat0: Failed loading <totcloudfrac>"
271  totcloudfrac(:)=0.5
272else
273  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
274             minval(totcloudfrac), maxval(totcloudfrac)
275endif
276
277! Height of oceanic ice (added by BC 2010)
278call get_field("hice",hice,found,indextime)
279if (.not.found) then
280  write(*,*) "phyetat0: Failed loading <hice>"
281!  call abort
282      do ig=1,ngrid
283      hice(ig)=0.
284      enddo
285else
286  write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
287             minval(hice), maxval(hice)
288endif
289
290! SLAB OCEAN (added by BC 2014)
291! nature of the surface
292call get_field("rnat",rnat,found,indextime)
293if (.not.found) then
294  write(*,*) "phyetat0: Failed loading <rnat>"
295      do ig=1,ngrid
296        rnat(ig)=1.
297      enddo
298else
299      do ig=1,ngrid
300        if((nint(rnat(ig)).eq.2).or.(nint(rnat(ig)).eq.0))then
301          rnat(ig)=0.
302        else
303          rnat(ig)=1.
304        endif     
305      enddo
306
307  write(*,*) "phyetat0: Nature of surface <rnat> range:", &
308             minval(rnat), maxval(rnat)
309endif
310! Pourcentage of sea ice cover
311call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
312if (.not.found) then
313  write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
314      do ig=1,ngrid
315      pctsrf_sic(ig)=0.
316      enddo
317else
318  write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
319             minval(pctsrf_sic), maxval(pctsrf_sic)
320endif
321! Slab ocean temperature (2 layers)
322call get_field("tslab",tslab,found,indextime)
323if (.not.found) then
324  write(*,*) "phyetat0: Failed loading <tslab>"
325      do ig=1,ngrid
326      do iq=1,noceanmx
327      tslab(ig,iq)=tsurf(ig)
328      enddo
329      enddo
330else
331  write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
332             minval(tslab), maxval(tslab)
333endif
334! Oceanic ice temperature
335call get_field("tsea_ice",tsea_ice,found,indextime)
336if (.not.found) then
337  write(*,*) "phyetat0: Failed loading <tsea_ice>"
338      do ig=1,ngrid
339      tsea_ice(ig)=273.15-1.8
340      enddo
341else
342  write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
343             minval(tsea_ice), maxval(tsea_ice)
344endif
345!  Oceanic ice quantity (kg/m^2)
346call get_field("sea_ice",sea_ice,found,indextime)
347if (.not.found) then
348  write(*,*) "phyetat0: Failed loading <sea_ice>"
349      do ig=1,ngrid
350      tsea_ice(ig)=0.
351      enddo
352else
353  write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
354             minval(sea_ice), maxval(sea_ice)
355endif
356
357
358
359
360! pbl wind variance
361call get_field("q2",q2,found,indextime)
362if (.not.found) then
363  write(*,*) "phyetat0: Failed loading <q2>"
364  q2(:,:)=0.001
365else
366  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
367             minval(q2), maxval(q2)
368endif
369
370! tracer on surface
371if (nq.ge.1) then
372  do iq=1,nq
373    txt=tname(iq)
374    if (txt.eq."h2o_vap") then
375      ! There is no surface tracer for h2o_vap;
376      ! "h2o_ice" should be loaded instead
377      txt="h2o_ice"
378      write(*,*) 'phyetat0: loading surface tracer', &
379                           ' h2o_ice instead of h2o_vap'
380    endif
381    call get_field(txt,qsurf(:,iq),found,indextime)
382    if (.not.found) then
383      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
384      write(*,*) "         ",trim(txt)," is set to zero"
385    else
386      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
387                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
388    endif
389  enddo
390endif ! of if (nq.ge.1)
391
392
393! Call to soil_settings, in order to read soil temperatures,
394! as well as thermal inertia and volumetric heat capacity
395
396!ym -> needed for saturn ?
397!ym call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
398
399
400!
401! close file:
402!
403call close_startphy
404
405END SUBROUTINE phyetat0_academic
Note: See TracBrowser for help on using the repository browser.