source: LMDZ6/branches/Ocean_skin/libf/phylmd/phyetat0.F90 @ 4115

Last change on this file since 4115 was 4024, checked in by lguez, 3 years ago

Only give initial value when needed

Only give initial value to delta_sst and delta_sal when needed:
when we want to send them to the ocean.

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