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
RevLine 
[232]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
[298]13  use mod_phys_lmdz_para, only : is_master
[232]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
[270]86      logical :: found,found_file
[232]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:
[270]101call open_startphy(fichnom,found_file)
[232]102
[270]103! Ehouarn, if file not found, then call tabfi with nid_start==0
104if (.not.found_file) then
[298]105  if (is_master) write(*,*) 'phyetat0_academic: call tabfi with nid_start=0'
[270]106  call tabfi (ngrid,0,Lmodif,tab0,day_ini,lmax,p_rad, &
[232]107                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
[270]108else
109  ! possibility to modify tab_cntrl in tabfi
[298]110  if (is_master) write(*,*)
111  if (is_master) write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
[270]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
[232]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:
[270]173if (found_file) then
174  call get_field("phisfi",phisfi,found)
175else
176  found=.false.
177endif
[232]178if (.not.found) then
[298]179  if (is_master) write(*,*) "phyetat0: Failed loading <phisfi>"
[232]180  phisfi(:)=0
181else
[298]182  if (is_master) write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
[232]183             minval(phisfi), maxval(phisfi)
184endif
185
186! Load bare ground albedo:
[270]187if (found_file) then
188  call get_field("albedodat",albedodat,found)
189else
190  found=.false.
191endif
[232]192if (.not.found) then
[298]193  if (is_master) write(*,*) "phyetat0: Failed loading <albedodat>"
[232]194  do ig=1,ngrid
195    albedodat(ig)=0.
196  enddo
197else
[298]198  if (is_master) write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
[232]199             minval(albedodat), maxval(albedodat)
200endif
201
202! ZMEA
[270]203if (found_file) then
204  call get_field("ZMEA",zmea,found)
205else
206  found=.false.
207endif
[232]208if (.not.found) then
209    zmea(:)=0.
210else
[298]211  if (is_master) write(*,*) "phyetat0: <ZMEA> range:", &
[232]212             minval(zmea), maxval(zmea)
213endif
214
215! ZSTD
[270]216if (found_file) then
217  call get_field("ZSTD",zstd,found)
218else
219  found=.false.
220endif
[232]221if (.not.found) then
[298]222  if (is_master) write(*,*) "phyetat0: Failed loading <ZSTD>"
[232]223   zstd(:)=0.
224else
[298]225  if (is_master) write(*,*) "phyetat0: <ZSTD> range:", &
[232]226             minval(zstd), maxval(zstd)
227endif
228
229! ZSIG
[270]230if (found_file) then
231  call get_field("ZSIG",zsig,found)
232else
233  found=.false.
234endif
[232]235if (.not.found) then
[298]236  if (is_master) write(*,*) "phyetat0: Failed loading <ZSIG>"
[232]237  zsig(:)=0.
238else
[298]239  if (is_master) write(*,*) "phyetat0: <ZSIG> range:", &
[232]240             minval(zsig), maxval(zsig)
241endif
242
243! ZGAM
[270]244if (found_file) then
245  call get_field("ZGAM",zgam,found)
246else
247  found=.false.
248endif
[232]249if (.not.found) then
[298]250  if (is_master) write(*,*) "phyetat0: Failed loading <ZGAM>"
[232]251  zgam(:)=0.
252else
[298]253  if (is_master) write(*,*) "phyetat0: <ZGAM> range:", &
[232]254             minval(zgam), maxval(zgam)
255endif
256
257! ZTHE
[270]258if (found_file) then
259  call get_field("ZTHE",zthe,found)
260else
261  found=.false.
262endif
[232]263if (.not.found) then
[298]264  if (is_master) write(*,*) "phyetat0: Failed loading <ZTHE>"
[232]265  zthe(:)=0.
266else
[298]267  if (is_master) write(*,*) "phyetat0: <ZTHE> range:", &
[232]268             minval(zthe), maxval(zthe)
269endif
270
271! Surface temperature :
[270]272if (found_file) then
273  call get_field("tsurf",tsurf,found,indextime)
274else
275  found=.false.
276endif
[232]277if (.not.found) then
[270]278  !mi initialising tsurf with pt(:,1)
279  !tsurf(:)=175.0
[232]280else
[298]281  if (is_master) write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
[232]282             minval(tsurf), maxval(tsurf)
283endif
284
285! Surface emissivity
[270]286if (found_file) then
287  call get_field("emis",emis,found,indextime)
288else
289  found=.false.
290endif
[232]291if (.not.found) then
[298]292  if (is_master) write(*,*) "phyetat0: Failed loading <emis>"
[232]293  emis(:)=0.5
294else
[298]295  if (is_master) write(*,*) "phyetat0: Surface emissivity <emis> range:", &
[232]296             minval(emis), maxval(emis)
297endif
298
299! Cloud fraction (added by BC 2010)
[270]300if (found_file) then
301  call get_field("cloudfrac",cloudfrac,found,indextime)
302else
303  found=.false.
304endif
[232]305if (.not.found) then
[298]306  if (is_master) write(*,*) "phyetat0: Failed loading <cloudfrac>"
[232]307  cloudfrac(:,:)=0.
308else
[298]309  if (is_master) write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
[232]310             minval(cloudfrac), maxval(cloudfrac)
311endif
312
313! Total cloud fraction (added by BC 2010)
[270]314if (found_file) then
315  call get_field("totcloudfrac",totcloudfrac,found,indextime)
316else
317  found=.false.
318endif
[232]319if (.not.found) then
[298]320  if (is_master) write(*,*) "phyetat0: Failed loading <totcloudfrac>"
[232]321  totcloudfrac(:)=0.5
322else
[298]323  if (is_master) write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
[232]324             minval(totcloudfrac), maxval(totcloudfrac)
325endif
326
327! Height of oceanic ice (added by BC 2010)
[270]328if (found_file) then
329  call get_field("hice",hice,found,indextime)
330else
331  found=.false.
332endif
[232]333if (.not.found) then
[298]334  if (is_master) write(*,*) "phyetat0: Failed loading <hice>"
[232]335!  call abort
336      do ig=1,ngrid
337      hice(ig)=0.
338      enddo
339else
[298]340  if (is_master) write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
[232]341             minval(hice), maxval(hice)
342endif
343
344! SLAB OCEAN (added by BC 2014)
345! nature of the surface
[270]346if (found_file) then
347  call get_field("rnat",rnat,found,indextime)
348else
349  found=.false.
350endif
[232]351if (.not.found) then
[298]352  if (is_master) write(*,*) "phyetat0: Failed loading <rnat>"
[232]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
[298]365  if (is_master) write(*,*) "phyetat0: Nature of surface <rnat> range:", &
[232]366             minval(rnat), maxval(rnat)
367endif
368! Pourcentage of sea ice cover
[270]369if (found_file) then
370  call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
371else
372  found=.false.
373endif
[232]374if (.not.found) then
[298]375  if (is_master) write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
[232]376      do ig=1,ngrid
377      pctsrf_sic(ig)=0.
378      enddo
379else
[298]380  if (is_master) write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
[232]381             minval(pctsrf_sic), maxval(pctsrf_sic)
382endif
383! Slab ocean temperature (2 layers)
[270]384if (found_file) then
385  call get_field("tslab",tslab,found,indextime)
386else
387  found=.false.
388endif
[232]389if (.not.found) then
[298]390  if (is_master) write(*,*) "phyetat0: Failed loading <tslab>"
[232]391      do ig=1,ngrid
392      do iq=1,noceanmx
393      tslab(ig,iq)=tsurf(ig)
394      enddo
395      enddo
396else
[298]397  if (is_master) write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
[232]398             minval(tslab), maxval(tslab)
399endif
400! Oceanic ice temperature
[270]401if (found_file) then
402  call get_field("tsea_ice",tsea_ice,found,indextime)
403else
404  found=.false.
405endif
[232]406if (.not.found) then
[298]407  if (is_master) write(*,*) "phyetat0: Failed loading <tsea_ice>"
[232]408      do ig=1,ngrid
409      tsea_ice(ig)=273.15-1.8
410      enddo
411else
[298]412  if (is_master) write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
[232]413             minval(tsea_ice), maxval(tsea_ice)
414endif
415!  Oceanic ice quantity (kg/m^2)
[270]416if (found_file) then
417  call get_field("sea_ice",sea_ice,found,indextime)
418else
419  found=.false.
420endif
[232]421if (.not.found) then
[298]422  if (is_master) write(*,*) "phyetat0: Failed loading <sea_ice>"
[232]423      do ig=1,ngrid
424      tsea_ice(ig)=0.
425      enddo
426else
[298]427  if (is_master) write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
[232]428             minval(sea_ice), maxval(sea_ice)
429endif
430
431
432
433
434! pbl wind variance
[270]435if (found_file) then
436  call get_field("q2",q2,found,indextime)
437else
438  found=.false.
439endif
[232]440if (.not.found) then
[298]441  if (is_master) write(*,*) "phyetat0: Failed loading <q2>"
[232]442  q2(:,:)=0.001
443else
[298]444  if (is_master) write(*,*) "phyetat0: PBL wind variance <q2> range:", &
[232]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"
[298]456      if (is_master) write(*,*) 'phyetat0: loading surface tracer', &
[232]457                           ' h2o_ice instead of h2o_vap'
458    endif
[270]459    if (found_file) then
460      call get_field(txt,qsurf(:,iq),found,indextime)
461    else
462      found=.false.
463    endif
[232]464    if (.not.found) then
[298]465      if (is_master) write(*,*) "phyetat0: Failed loading <",trim(txt),">"
466      if (is_master) write(*,*) "         ",trim(txt)," is set to zero"
[232]467    else
[298]468      if (is_master) write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
[232]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.