source: LMDZ5/trunk/libf/phylmd/phyetat0.F90 @ 2989

Last change on this file since 2989 was 2979, checked in by Ehouarn Millour, 7 years ago

Missing initializations that caused problems when running in aquaplanet mode.
EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.9 KB
Line 
1! $Id: phyetat0.F90 2979 2017-08-03 08:25:17Z musat $
2
3SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
4
5  USE dimphy, only: klon, zmasq, klev
6  USE iophy, ONLY : init_iophy_new
7  USE ocean_cpl_mod,    ONLY : ocean_cpl_init
8  USE fonte_neige_mod,  ONLY : fonte_neige_init
9  USE pbl_surface_mod,  ONLY : pbl_surface_init
10  USE surface_data,     ONLY : type_ocean, version_ocean
11  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
12       qsol, fevap, z0m, z0h, agesno, &
13       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
14       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, &
15       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, radpas, radsol, rain_fall, ratqs, &
16       rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
17       solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
18       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
19       wake_s, wake_dens, zgam, zmax0, zmea, zpic, zsig, &
20       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg
21!FC
22  USE geometry_mod, ONLY : longitude_deg, latitude_deg
23  USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
24  USE infotrac_phy, only: nbtr, nqo, type_trac, tname, niadv
25  USE traclmdz_mod,    ONLY : traclmdz_from_restart
26  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
27  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
28  USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init
29  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
30
31  IMPLICIT none
32  !======================================================================
33  ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
34  ! Objet: Lecture de l'etat initial pour la physique
35  !======================================================================
36  include "netcdf.inc"
37  include "dimsoil.h"
38  include "clesphys.h"
39  include "thermcell.h"
40  include "compbl.h"
41  include "YOMCST.h"
42  !======================================================================
43  CHARACTER*(*) fichnom
44
45  ! les variables globales lues dans le fichier restart
46
47  REAL tsoil(klon, nsoilmx, nbsrf)
48  REAL qsurf(klon, nbsrf)
49  REAL snow(klon, nbsrf)
50  real fder(klon)
51  REAL run_off_lic_0(klon)
52  REAL fractint(klon)
53  REAL trs(klon, nbtr)
54  REAL zts(klon)
55  ! pour drag arbres FC
56  REAL drg_ter(klon,klev)
57
58  CHARACTER*6 ocean_in
59  LOGICAL ok_veget_in
60
61  INTEGER        longcles
62  PARAMETER    ( longcles = 20 )
63  REAL clesphy0( longcles )
64
65  REAL xmin, xmax
66
67  INTEGER nid, nvarid
68  INTEGER ierr, i, nsrf, isoil , k
69  INTEGER length
70  PARAMETER (length=100)
71  INTEGER it, iiq, isw
72  REAL tab_cntrl(length), tabcntr0(length)
73  CHARACTER*7 str7
74  CHARACTER*2 str2
75  LOGICAL :: found,phyetat0_get,phyetat0_srf
76  REAL :: lon_startphy(klon), lat_startphy(klon)
77
78  ! FH1D
79  !     real iolat(jjm+1)
80  !real iolat(jjm+1-1/(iim*jjm))
81
82  ! Ouvrir le fichier contenant l'etat initial:
83
84  CALL open_startphy(fichnom)
85
86  ! Lecture des parametres de controle:
87
88  CALL get_var("controle", tab_cntrl)
89
90!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
91  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
92  ! Les constantes de la physiques sont lues dans la physique seulement.
93  ! Les egalites du type
94  !             tab_cntrl( 5 )=clesphy0(1)
95  ! sont remplacees par
96  !             clesphy0(1)=tab_cntrl( 5 )
97  ! On inverse aussi la logique.
98  ! On remplit les tab_cntrl avec les parametres lus dans les .def
99!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
100
101  DO i = 1, length
102     tabcntr0( i ) = tab_cntrl( i )
103  ENDDO
104
105  tab_cntrl(1)=pdtphys
106  tab_cntrl(2)=radpas
107
108  ! co2_ppm : value from the previous time step
109  IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
110     co2_ppm = tab_cntrl(3)
111     RCO2    = co2_ppm * 1.0e-06  * 44.011/28.97
112     ! ELSE : keep value from .def
113  END IF
114
115  ! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
116  co2_ppm0   = tab_cntrl(16)
117
118  solaire_etat0      = tab_cntrl(4)
119  tab_cntrl(5)=iflag_con
120  tab_cntrl(6)=nbapp_rad
121
122  if (cycle_diurne) tab_cntrl( 7) =1.
123  if (soil_model) tab_cntrl( 8) =1.
124  if (new_oliq) tab_cntrl( 9) =1.
125  if (ok_orodr) tab_cntrl(10) =1.
126  if (ok_orolf) tab_cntrl(11) =1.
127  if (ok_limitvrai) tab_cntrl(12) =1.
128
129  itau_phy = tab_cntrl(15)
130
131  clesphy0(1)=tab_cntrl( 5 )
132  clesphy0(2)=tab_cntrl( 6 )
133  clesphy0(3)=tab_cntrl( 7 )
134  clesphy0(4)=tab_cntrl( 8 )
135  clesphy0(5)=tab_cntrl( 9 )
136  clesphy0(6)=tab_cntrl( 10 )
137  clesphy0(7)=tab_cntrl( 11 )
138  clesphy0(8)=tab_cntrl( 12 )
139
140  ! set time iteration
141   CALL init_iteration(itau_phy)
142
143  ! read latitudes and make a sanity check (because already known from dyn)
144  CALL get_field("latitude",lat_startphy)
145  DO i=1,klon
146    IF (ABS(lat_startphy(i)-latitude_deg(i))>=1) THEN
147      WRITE(*,*) "phyetat0: Error! Latitude discrepancy wrt startphy file:",&
148                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
149                 " latitude_deg(i)=",latitude_deg(i)
150      ! This is presumably serious enough to abort run
151      CALL abort_physic("phyetat0","discrepancy in latitudes!",1)
152    ENDIF
153    IF (ABS(lat_startphy(i)-latitude_deg(i))>=0.0001) THEN
154      WRITE(*,*) "phyetat0: Warning! Latitude discrepancy wrt startphy file:",&
155                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
156                 " latitude_deg(i)=",latitude_deg(i)
157    ENDIF
158  ENDDO
159
160  ! read longitudes and make a sanity check (because already known from dyn)
161  CALL get_field("longitude",lon_startphy)
162  DO i=1,klon
163    IF (ABS(lon_startphy(i)-longitude_deg(i))>=1) THEN
164      WRITE(*,*) "phyetat0: Error! Longitude discrepancy wrt startphy file:",&
165                 " i=",i," lon_startphy(i)=",lon_startphy(i),&
166                 " longitude_deg(i)=",longitude_deg(i)
167      ! This is presumably serious enough to abort run
168      CALL abort_physic("phyetat0","discrepancy in longitudes!",1)
169    ENDIF
170    IF (ABS(lon_startphy(i)-longitude_deg(i))>=0.0001) THEN
171      WRITE(*,*) "phyetat0: Warning! Longitude discrepancy wrt startphy file:",&
172                 " i=",i," lon_startphy(i)=",lon_startphy(i),&
173                 " longitude_deg(i)=",longitude_deg(i)
174    ENDIF
175  ENDDO
176
177  ! Lecture du masque terre mer
178
179  CALL get_field("masque", zmasq, found)
180  IF (.NOT. found) THEN
181     PRINT*, 'phyetat0: Le champ <masque> est absent'
182     PRINT *, 'fichier startphy non compatible avec phyetat0'
183  ENDIF
184
185  ! Lecture des fractions pour chaque sous-surface
186
187  ! initialisation des sous-surfaces
188
189  pctsrf = 0.
190
191  ! fraction de terre
192
193  CALL get_field("FTER", pctsrf(:, is_ter), found)
194  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FTER> est absent'
195
196  ! fraction de glace de terre
197
198  CALL get_field("FLIC", pctsrf(:, is_lic), found)
199  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FLIC> est absent'
200
201  ! fraction d'ocean
202
203  CALL get_field("FOCE", pctsrf(:, is_oce), found)
204  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FOCE> est absent'
205
206  ! fraction glace de mer
207
208  CALL get_field("FSIC", pctsrf(:, is_sic), found)
209  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FSIC> est absent'
210
211  !  Verification de l'adequation entre le masque et les sous-surfaces
212
213  fractint( 1 : klon) = pctsrf(1 : klon, is_ter)  &
214       + pctsrf(1 : klon, is_lic)
215  DO i = 1 , klon
216     IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
217        WRITE(*, *) 'phyetat0: attention fraction terre pas ',  &
218             'coherente ', i, zmasq(i), pctsrf(i, is_ter) &
219             , pctsrf(i, is_lic)
220        WRITE(*, *) 'Je force la coherence zmasq=fractint'
221        zmasq(i) = fractint(i)
222     ENDIF
223  END DO
224  fractint (1 : klon) =  pctsrf(1 : klon, is_oce)  &
225       + pctsrf(1 : klon, is_sic)
226  DO i = 1 , klon
227     IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
228        WRITE(*, *) 'phyetat0 attention fraction ocean pas ',  &
229             'coherente ', i, zmasq(i) , pctsrf(i, is_oce) &
230             , pctsrf(i, is_sic)
231        WRITE(*, *) 'Je force la coherence zmasq=1.-fractint'
232        zmasq(i) = 1. - fractint(i)
233     ENDIF
234  END DO
235
236!===================================================================
237! Lecture des temperatures du sol:
238!===================================================================
239
240  found=phyetat0_get(1,ftsol(:,1),"TS","Surface temperature",283.)
241  IF (found) THEN
242     DO nsrf=2,nbsrf
243        ftsol(:,nsrf)=ftsol(:,1)
244     ENDDO
245  ELSE
246     found=phyetat0_srf(1,ftsol,"TS","Surface temperature",283.)
247  ENDIF
248
249!===================================================================
250  ! Lecture des albedo difus et direct
251!===================================================================
252
253  DO nsrf = 1, nbsrf
254     DO isw=1, nsw
255        IF (isw.GT.99) THEN
256           PRINT*, "Trop de bandes SW"
257           call abort_physic("phyetat0", "", 1)
258        ENDIF
259        WRITE(str2, '(i2.2)') isw
260        found=phyetat0_srf(1,falb_dir(:, isw,:),"A_dir_SW"//str2//"srf","Direct Albedo",0.2)
261        found=phyetat0_srf(1,falb_dif(:, isw,:),"A_dif_SW"//str2//"srf","Direct Albedo",0.2)
262     ENDDO
263  ENDDO
264
265  found=phyetat0_srf(1,u10m,"U10M","u a 10m",0.)
266  found=phyetat0_srf(1,v10m,"V10M","v a 10m",0.)
267
268!===================================================================
269  ! Lecture des temperatures du sol profond:
270!===================================================================
271
272   DO isoil=1, nsoilmx
273        IF (isoil.GT.99) THEN
274           PRINT*, "Trop de couches "
275           call abort_physic("phyetat0", "", 1)
276        ENDIF
277        WRITE(str2,'(i2.2)') isoil
278        found=phyetat0_srf(1,tsoil(:, isoil,:),"Tsoil"//str2//"srf","Temp soil",0.)
279        IF (.NOT. found) THEN
280           PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
281           PRINT*, "          Il prend donc la valeur de surface"
282           tsoil(:, isoil, :)=ftsol(:, :)
283        ENDIF
284   ENDDO
285
286!=======================================================================
287! Lecture precipitation/evaporation
288!=======================================================================
289
290  found=phyetat0_srf(1,qsurf,"QS","Near surface hmidity",0.)
291  found=phyetat0_get(1,qsol,"QSOL","Surface hmidity / bucket",0.)
292  found=phyetat0_srf(1,snow,"SNOW","Surface snow",0.)
293  found=phyetat0_srf(1,fevap,"EVAP","evaporation",0.)
294  found=phyetat0_get(1,snow_fall,"snow_f","snow fall",0.)
295  found=phyetat0_get(1,rain_fall,"rain_f","rain fall",0.)
296
297!=======================================================================
298! Radiation
299!=======================================================================
300
301  found=phyetat0_get(1,solsw,"solsw","net SW radiation surf",0.)
302  found=phyetat0_get(1,sollw,"sollw","net LW radiation surf",0.)
303  found=phyetat0_get(1,sollwdown,"sollwdown","down LW radiation surf",0.)
304  IF (.NOT. found) THEN
305     sollwdown = 0. ;  zts=0.
306     do nsrf=1,nbsrf
307        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
308     enddo
309     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
310  ENDIF
311
312  found=phyetat0_get(1,radsol,"RADS","Solar radiation",0.)
313  found=phyetat0_get(1,fder,"fder","Flux derivative",0.)
314
315
316  ! Lecture de la longueur de rugosite
317  found=phyetat0_srf(1,z0m,"RUG","Z0m ancien",0.001)
318  IF (found) THEN
319     z0h(:,1:nbsrf)=z0m(:,1:nbsrf)
320  ELSE
321     found=phyetat0_srf(1,z0m,"Z0m","Roughness length, momentum ",0.001)
322     found=phyetat0_srf(1,z0h,"Z0h","Roughness length, enthalpy ",0.001)
323  ENDIF
324!FC
325  IF (ifl_pbltree>0) THEN
326!CALL get_field("FTER", pctsrf(:, is_ter), found)
327    treedrg(:,1:klev,1:nbsrf)= 0.0
328    CALL get_field("treedrg_ter", drg_ter(:,:), found)
329!  found=phyetat0_srf(1,treedrg,"treedrg","drag from vegetation" , 0.)
330    !lecture du profile de freinage des arbres
331    IF (.not. found ) THEN
332      treedrg(:,1:klev,1:nbsrf)= 0.0
333    ELSE
334      treedrg(:,1:klev,is_ter)= drg_ter(:,:)
335!     found=phyetat0_srf(klev,treedrg,"treedrg","freinage arbres",0.)
336    ENDIF
337  ELSE
338    ! initialize treedrg(), because it will be written in restartphy.nc
339    treedrg(:,:,:) = 0.0
340  ENDIF
341
342  ! Lecture de l'age de la neige:
343  found=phyetat0_srf(1,agesno,"AGESNO","SNOW AGE",0.001)
344
345  ancien_ok=.true.
346  ancien_ok=ancien_ok.AND.phyetat0_get(klev,t_ancien,"TANCIEN","TANCIEN",0.)
347  ancien_ok=ancien_ok.AND.phyetat0_get(klev,q_ancien,"QANCIEN","QANCIEN",0.)
348  ancien_ok=ancien_ok.AND.phyetat0_get(klev,ql_ancien,"QLANCIEN","QLANCIEN",0.)
349  ancien_ok=ancien_ok.AND.phyetat0_get(klev,qs_ancien,"QSANCIEN","QSANCIEN",0.)
350  ancien_ok=ancien_ok.AND.phyetat0_get(klev,u_ancien,"UANCIEN","UANCIEN",0.)
351  ancien_ok=ancien_ok.AND.phyetat0_get(klev,v_ancien,"VANCIEN","VANCIEN",0.)
352  ancien_ok=ancien_ok.AND.phyetat0_get(1,prw_ancien,"PRWANCIEN","PRWANCIEN",0.)
353  ancien_ok=ancien_ok.AND.phyetat0_get(1,prlw_ancien,"PRLWANCIEN","PRLWANCIEN",0.)
354  ancien_ok=ancien_ok.AND.phyetat0_get(1,prsw_ancien,"PRSWANCIEN","PRSWANCIEN",0.)
355
356  ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
357  !          dummy values (as is the case when generated by ce0l,
358  !          or by iniaqua)
359  if ( (maxval(q_ancien).eq.minval(q_ancien))       .or. &
360       (maxval(ql_ancien).eq.minval(ql_ancien))     .or. &
361       (maxval(qs_ancien).eq.minval(qs_ancien))     .or. &
362       (maxval(prw_ancien).eq.minval(prw_ancien))   .or. &
363       (maxval(prlw_ancien).eq.minval(prlw_ancien)) .or. &
364       (maxval(prsw_ancien).eq.minval(prsw_ancien)) .or. &
365       (maxval(t_ancien).eq.minval(t_ancien)) ) then
366    ancien_ok=.false.
367  endif
368
369  found=phyetat0_get(klev,clwcon,"CLWCON","CLWCON",0.)
370  found=phyetat0_get(klev,rnebcon,"RNEBCON","RNEBCON",0.)
371  found=phyetat0_get(klev,ratqs,"RATQS","RATQS",0.)
372
373  found=phyetat0_get(1,run_off_lic_0,"RUNOFFLIC0","RUNOFFLIC0",0.)
374
375!==================================
376!  TKE
377!==================================
378!
379  IF (iflag_pbl>1) then
380     found=phyetat0_srf(klev+1,pbl_tke,"TKE","Turb. Kinetic. Energ. ",1.e-8)
381  ENDIF
382
383  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1 ) then
384    found=phyetat0_srf(klev+1,wake_delta_pbl_tke,"DELTATKE","Del TKE wk/env",0.)
385    found=phyetat0_srf(1,delta_tsurf,"DELTA_TSURF","Delta Ts wk/env ",0.)
386  ENDIF   !(iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1 )
387
388!==================================
389!  thermiques, poches, convection
390!==================================
391
392! Emanuel
393  found=phyetat0_get(klev,sig1,"sig1","sig1",0.)
394  found=phyetat0_get(klev,w01,"w01","w01",0.)
395
396! Wake
397  found=phyetat0_get(klev,wake_deltat,"WAKE_DELTAT","Delta T wake/env",0.)
398  found=phyetat0_get(klev,wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
399  found=phyetat0_get(1,wake_s,"WAKE_S","Wake frac. area",0.)
400  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",0.)
401  found=phyetat0_get(1,wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
402  found=phyetat0_get(1,wake_pe,"WAKE_PE","WAKE_PE",0.)
403  found=phyetat0_get(1,wake_fip,"WAKE_FIP","WAKE_FIP",0.)
404
405! Thermiques
406  found=phyetat0_get(1,zmax0,"ZMAX0","ZMAX0",40.)
407  found=phyetat0_get(1,f0,"F0","F0",1.e-5)
408  found=phyetat0_get(klev+1,fm_therm,"FM_THERM","Thermals mass flux",0.)
409  found=phyetat0_get(klev,entr_therm,"ENTR_THERM","Thermals Entrain.",0.)
410  found=phyetat0_get(klev,detr_therm,"DETR_THERM","Thermals Detrain.",0.)
411
412! ALE/ALP
413  found=phyetat0_get(1,ale_bl,"ALE_BL","ALE BL",0.)
414  found=phyetat0_get(1,ale_bl_trig,"ALE_BL_TRIG","ALE BL_TRIG",0.)
415  found=phyetat0_get(1,alp_bl,"ALP_BL","ALP BL",0.)
416
417!===========================================
418  ! Read and send field trs to traclmdz
419!===========================================
420
421  IF (type_trac == 'lmdz') THEN
422     DO it=1, nbtr                                                                 
423!!        iiq=niadv(it+2)                                                           ! jyg
424        iiq=niadv(it+nqo)                                                           ! jyg
425        found=phyetat0_get(1,trs(:,it),"trs_"//tname(iiq), &
426              "Surf trac"//tname(iiq),0.)
427     END DO
428     CALL traclmdz_from_restart(trs)
429
430     IF (carbon_cycle_cpl) THEN
431        ALLOCATE(co2_send(klon), stat=ierr)
432        IF (ierr /= 0) CALL abort_physic('phyetat0', 'pb allocation co2_send', 1)
433        found=phyetat0_get(1,co2_send,"co2_send","co2 send",0.)
434     END IF
435  END IF
436
437!===========================================
438!  ondes de gravite / relief
439!===========================================
440
441!  ondes de gravite non orographiques
442  if (ok_gwd_rando) found = &
443       phyetat0_get(klev,du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.)
444  IF (.not. ok_hines .and. ok_gwd_rando) found &
445       = phyetat0_get(klev,du_gwd_front,"du_gwd_front","du_gwd_front",0.)
446
447!  prise en compte du relief sous-maille
448  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
449  found=phyetat0_get(1,zstd,"ZSTD","sub grid orography",0.)
450  found=phyetat0_get(1,zsig,"ZSIG","sub grid orography",0.)
451  found=phyetat0_get(1,zgam,"ZGAM","sub grid orography",0.)
452  found=phyetat0_get(1,zthe,"ZTHE","sub grid orography",0.)
453  found=phyetat0_get(1,zpic,"ZPIC","sub grid orography",0.)
454  found=phyetat0_get(1,zval,"ZVAL","sub grid orography",0.)
455  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
456  found=phyetat0_get(1,rugoro,"RUGSREL","sub grid orography",0.)
457
458!===========================================
459! Initialize ocean
460!===========================================
461
462  IF ( type_ocean == 'slab' ) THEN
463      CALL ocean_slab_init(dtime, pctsrf)
464      IF (nslay.EQ.1) THEN
465        found=phyetat0_get(1,tslab,"tslab01","tslab",0.)
466        IF (.NOT. found) THEN
467            found=phyetat0_get(1,tslab,"tslab","tslab",0.)
468        END IF
469      ELSE
470          DO i=1,nslay
471            WRITE(str2,'(i2.2)') i
472            found=phyetat0_get(1,tslab(:,i),"tslab"//str2,"tslab",0.) 
473          END DO
474      END IF
475      IF (.NOT. found) THEN
476          PRINT*, "phyetat0: Le champ <tslab> est absent"
477          PRINT*, "Initialisation a tsol_oce"
478          DO i=1,nslay
479              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
480          END DO
481      END IF
482
483      ! Sea ice variables
484      IF (version_ocean == 'sicINT') THEN
485          found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
486          IF (.NOT. found) THEN
487              PRINT*, "phyetat0: Le champ <tice> est absent"
488              PRINT*, "Initialisation a tsol_sic"
489                  tice(:)=ftsol(:,is_sic)
490          END IF
491          found=phyetat0_get(1,seaice,"seaice","seaice",0.)
492          IF (.NOT. found) THEN
493              PRINT*, "phyetat0: Le champ <seaice> est absent"
494              PRINT*, "Initialisation a 0/1m suivant fraction glace"
495              seaice(:)=0.
496              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
497                  seaice=917.
498              END WHERE
499          END IF
500      END IF !sea ice INT
501  END IF ! Slab       
502
503  ! on ferme le fichier
504  CALL close_startphy
505
506  ! Initialize module pbl_surface_mod
507
508  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
509
510  ! Initialize module ocean_cpl_mod for the case of coupled ocean
511  IF ( type_ocean == 'couple' ) THEN
512     CALL ocean_cpl_init(dtime, longitude_deg, latitude_deg)
513  ENDIF
514
515  CALL init_iophy_new(latitude_deg, longitude_deg)
516
517  ! Initilialize module fonte_neige_mod     
518  CALL fonte_neige_init(run_off_lic_0)
519
520END SUBROUTINE phyetat0
521
522!===================================================================
523FUNCTION phyetat0_get(nlev,field,name,descr,default)
524!===================================================================
525! Lecture d'un champ avec contrôle
526! Function logique dont le resultat indique si la lecture
527! s'est bien passée
528! On donne une valeur par defaut dans le cas contraire
529!===================================================================
530
531USE iostart, ONLY : get_field
532USE dimphy, only: klon
533USE print_control_mod, ONLY: lunout
534
535IMPLICIT NONE
536
537LOGICAL phyetat0_get
538
539! arguments
540INTEGER,INTENT(IN) :: nlev
541CHARACTER*(*),INTENT(IN) :: name,descr
542REAL,INTENT(IN) :: default
543REAL,DIMENSION(klon,nlev),INTENT(INOUT) :: field
544
545! Local variables
546LOGICAL found
547
548   CALL get_field(name, field, found)
549   IF (.NOT. found) THEN
550     WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
551     WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
552     field(:,:)=default
553   ENDIF
554   WRITE(lunout,*) name, descr, MINval(field),MAXval(field)
555   phyetat0_get=found
556
557RETURN
558END FUNCTION phyetat0_get
559
560!================================================================
561FUNCTION phyetat0_srf(nlev,field,name,descr,default)
562!===================================================================
563! Lecture d'un champ par sous-surface avec contrôle
564! Function logique dont le resultat indique si la lecture
565! s'est bien passée
566! On donne une valeur par defaut dans le cas contraire
567!===================================================================
568
569USE iostart, ONLY : get_field
570USE dimphy, only: klon
571USE indice_sol_mod, only: nbsrf
572USE print_control_mod, ONLY: lunout
573
574IMPLICIT NONE
575
576LOGICAL phyetat0_srf
577! arguments
578INTEGER,INTENT(IN) :: nlev
579CHARACTER*(*),INTENT(IN) :: name,descr
580REAL,INTENT(IN) :: default
581REAL,DIMENSION(klon,nlev,nbsrf),INTENT(INOUT) :: field
582
583! Local variables
584LOGICAL found,phyetat0_get
585INTEGER nsrf
586CHARACTER*2 str2
587 
588     IF (nbsrf.GT.99) THEN
589        WRITE(lunout,*) "Trop de sous-mailles"
590        call abort_physic("phyetat0", "", 1)
591     ENDIF
592
593     DO nsrf = 1, nbsrf
594        WRITE(str2, '(i2.2)') nsrf
595        found= phyetat0_get(nlev,field(:,:, nsrf), &
596        name//str2,descr//" srf:"//str2,default)
597     ENDDO
598
599     phyetat0_srf=found
600
601RETURN
602END FUNCTION phyetat0_srf
603
Note: See TracBrowser for help on using the repository browser.