source: LMDZ6/trunk/libf/phylmd/phyetat0.F90 @ 3513

Last change on this file since 3513 was 3505, checked in by yann meurdesoif, 6 years ago

Solve some wrong discrepency problem when comparing longitude from a restartphy file. The current discrepency test is not detecting that 360°==0°, so in some case it may stop the run for a wrong reason.

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