source: LMDZ6/trunk/libf/phylmdiso/phyetat0.F90 @ 3932

Last change on this file since 3932 was 3927, checked in by Laurent Fairhead, 3 years ago

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