source: LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90 @ 5278

Last change on this file since 5278 was 5274, checked in by abarral, 23 hours ago

Replace yomcst.h by existing module

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