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

Last change on this file since 1303 was 1297, checked in by bclmd, 10 years ago

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