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

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

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