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

Last change on this file since 3458 was 3458, checked in by lguez, 6 years ago

Introduce variable activate_ocean_skin in module config_ocean_skin_m.

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