source: LMDZ6/trunk/libf/phylmd/sisvat/surf_sisvat_mod.F90 @ 3761

Last change on this file since 3761 was 3518, checked in by Laurent Fairhead, 6 years ago

Modifications needed as put_field call has changed

  • 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
File size: 84.7 KB
RevLine 
[1865]1MODULE surf_sisvat_mod
2  USE dimphy
3  IMPLICIT NONE
4
5  INTEGER, PARAMETER    :: nsnowmx=35
6  INTEGER, PARAMETER    :: nsismx=46         ! = nsnowmx + nsoilmx
7
8CONTAINS
9
10                                       
11
12  SUBROUTINE surf_sisvat(knon,rlon,rlat, ikl2i, itime, dtime, debut, lafin, &
13             rmu0, swdown, lwdown, pexner, ps, p1lay, &
14             precip_rain, precip_snow, precip_snow_adv, snow_adv, &
15             bl_height, wind_velo, temp_air, dens_air, spechum, tsurf, &
16             rugos, snow_cont_air, alb_soil, slope, cloudf, &
17             radsol, qsol, tsoil, snow, snowhgt, qsnow, to_ice, sissnow, agesno, &
18             AcoefH, AcoefQ, BcoefH, BcoefQ, cdragh, &
19             runoff_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &     
20             tsurf_new, alb1, alb2, alb3, &
21             emis_new, z0_new, qsurf)       
22                                                                             
23! +------------------------------------------------------------------------+   
24! |                                                                        |   
25! |   SubRoutine surf_sisvat: Interface between LMDZ and landice scheme    |
26! |     of the SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)|   
27! |                                                                        |
28! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
29! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
30! |                                                                        |   
31! +------------------------------------------------------------------------+   
32! |   
33! |   In the current setup, SISVAT is used only to model the land ice      |
34! |   part of the surface; hence it is called with the compressed variables|
35! |   from pbl_surface, and only by the surf_landice routine.              |
36! |                                                                        |   
37! |   In this interface it is assumed that the partitioning of the soil,   |
38! |   and hence the number of grid points is constant during a simulation, |
39! |   hence eg. snow properties remain stored in the global SISVAT         |
40! |   variables between the calls and don't need to be handed over as      |
41! |   arguments. When the partitioning is supposed to change, make sure to |
42! |   update the variables.                                                |
43! |                                                                        | 
44! |   INPUT                                                                | 
45! |   ^^^^^     VegMod: SISVAT    is set up when .T.                       | 
46! |             SnoMod: Snow Pack is set up when .T.                       | 
47! |             reaLBC: Update Bound.Condit.when .T.                       |   
48! |                                                                        | 
49! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV)                        | 
50! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |   
51! |                                                                        | 
52! |   Preprocessing  Option: SISVAT PHYSICS                                |   
53! |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^^^^^^                                |   
54! | #                       #HY                                            |   
55! | #                       #SN: Snow         Model                        |   
56! | #                       #BS: Blowing Snow Parameterization             |   
57! |                                                                        |   
58! | #                       #DS: diffuse radiation differing from direct   |   
59! |                             (variable RADsod must still be included)   |
60! | #                       #CP: SBL,                       Col de Porte   |   
61! | #                       #cp  Solar Radiation,           Col de Porte   |   
62! | #                       #AG: Snow Ageing,               Col de Porte   |   
63! |                                                                        |   
64! |   Preprocessing  Option: SISVAT IO                                     |   
65! |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^                                     |   
66! |   FILE                 |      CONTENT                                  |   
67! |   ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |   
68! | # ANI.yyyymmdd.LAB.nc  | #NC: OUTPUT on NetCDF File (Stand Alone EXP.) |   
69! | # SISVAT_iii_jjj_n     | #WV: OUTPUT on ASCII  File (SISVAT Variables) |   
70! | # SISVATtroubles       | #VF: OUTPUT on ASCII  File (SISVAT Troubles)  |   
71! |                        |                                               |   
72! | #                      | #ES: OUTPUT/Verification: Energy Conservation |   
73! | #                      | #E2: OUTPUT/Verification: Energy Consrv.2e pt.|   
74! |                        |                           (no premature stop) |   
75! | #                      | #MW: OUTPUT/Verification: H2O    Conservation |   
76! | #                      | #MS: OUTPUT/Verification: * Mass Conservation |   
77! | #                      | #MI: OUTPUT/Verification: SeaIce Conservation |   
78! |                        |                                               |   
79! | # SISVAT__zSn.OUT      | #as: OUTPUT/Verification: Snow Layers Agrega. |   
80! | # SISVAT__SnO.OUT      | #aw: OUTPUT/Verification: Albedo Parameteriz. |   
81! | # SISVATe_qSn.OUT      | #em: OUTPUT/Verification: Energy/Water Budget |   
82! | # SISVATu_qSn.OUT      | #su: OUTPUT/Verification: Slush  Parameteriz. |   
83! | # SISVATw_qSo.OUT      | #mw: OUTPUT/Verif+Detail: H2O    Conservation |   
84! |                        |                                               |   
85! | # SISVAT__GSn.OUT      | #VP: OUTPUT/Verification: Snow   Properties   |   
86! | # SISVAT__wEq.OUT      | #EQ: OUTPUT/Verification: Snow/Ice Water Eqv. |   
87! |                        |                                               |   
88! | #                      | #VR: VERIFICATION OUTPUT                      |   
89! | #                      | #WR: Additional   OUTPUT                      |   
90! |                                                                        |   
91! +------------------------------------------------------------------------+
92                                                       
93    USE VAR_SV
94    USE VARdSV
95    USE VARxSV
96    USE VARySV
97    USE VARtSV
98 
99    USE VARdCP
100    USE VARphy, ra_earth=>ra
[1990]101    USE YOMCST_SISVAT
[1865]102                                       
103    IMPLICIT NONE   
104                                                               
105! +--INTERFACE Variables                                                     
106! +  ===================
107                                           
108    include  "dimsoil.h"                                     
109                                                                   
110
111! +--Global Variables                                                         
112! +  ================ 
113! Input Variables for SISVAT
114    INTEGER,               INTENT(IN)      :: knon
115    INTEGER,               INTENT(IN)      :: itime   
116    REAL,                  INTENT(IN)      :: dtime
117    LOGICAL,               INTENT(IN)      :: debut     ! true if first step
118    LOGICAL,               INTENT(IN)      :: lafin     ! true if last step
119
120    INTEGER, DIMENSION(klon), INTENT(IN)   :: ikl2i     ! Index Decompression
121    REAL, DIMENSION(klon), INTENT(IN)      :: rlon, rlat
122    REAL, DIMENSION(klon), INTENT(IN)      :: rmu0      ! cos sol. zenith angle
123    REAL, DIMENSION(klon), INTENT(IN)      :: swdown    !
124    REAL, DIMENSION(klon), INTENT(IN)      :: lwdown    !
125    REAL, DIMENSION(klon), INTENT(IN)      :: pexner    ! Exner potential
126    REAL, DIMENSION(klon), INTENT(IN)      :: precip_rain, precip_snow
127    REAL, DIMENSION(klon), INTENT(IN)      :: precip_snow_adv, snow_adv
128                                                        !Snow Drift
129    REAL, DIMENSION(klon), INTENT(IN)      :: bl_height, wind_velo
130    REAL, DIMENSION(klon), INTENT(IN)      :: temp_air, spechum, ps,p1lay
131    REAL, DIMENSION(klon), INTENT(IN)      :: dens_air, tsurf           
132    REAL, DIMENSION(klon), INTENT(IN)      :: rugos,snow_cont_air
133    REAL, DIMENSION(klon), INTENT(IN)      :: alb_soil, slope
134    REAL, DIMENSION(klon), INTENT(IN)      :: cloudf   
135    REAL, DIMENSION(klon), INTENT(IN)      :: AcoefH, AcoefQ
136    REAL, DIMENSION(klon), INTENT(IN)      :: BcoefH, BcoefQ
137    REAL, DIMENSION(klon), INTENT(IN)      :: cdragh
138
139! Variables exchanged between LMDZ and SISVAT
140    REAL, DIMENSION(klon,nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
141    REAL, DIMENSION(klon), INTENT(OUT)     :: qsol      ! Soil Water Content   
142    REAL, DIMENSION(klon), INTENT(INOUT)   :: snow      ! Tot snow mass [kg/m2]
143    REAL, DIMENSION(klon), INTENT(IN)      :: radsol    ! Surface absorbed rad.
144
145! Output Variables from SISVAT
146    REAL, DIMENSION(klon), INTENT(OUT)     :: alb1      ! Albedo SW
147    REAL, DIMENSION(klon), INTENT(OUT)     :: alb2,alb3      ! Albedo LW
148    REAL, DIMENSION(klon), INTENT(OUT)     :: emis_new  ! Surface Emissivity
149    REAL, DIMENSION(klon), INTENT(OUT)     :: z0_new    ! Momentum Roughn Lgt
150    REAL, DIMENSION(klon), INTENT(OUT)     :: runoff_lic ! Runoff           
151    REAL, DIMENSION(klon), INTENT(OUT)     :: dflux_s   ! d/dT sens. ht flux 
152    REAL, DIMENSION(klon), INTENT(OUT)     :: dflux_l   ! d/dT latent ht flux
153    REAL, DIMENSION(klon), INTENT(OUT)     :: fluxsens  ! Sensible ht flux   
154    REAL, DIMENSION(klon), INTENT(OUT)     :: fluxlat   ! Latent heat flux
155    REAL, DIMENSION(klon), INTENT(OUT)     :: evap      ! Evaporation
156    REAL, DIMENSION(klon), INTENT(OUT)     :: agesno    ! Snow age (top layer)
157    REAL, DIMENSION(klon), INTENT(OUT)     :: tsurf_new ! Surface Temperature
158    REAL, DIMENSION(klon), INTENT(OUT)     :: qsurf     ! Surface Humidity
159    REAL, DIMENSION(klon), INTENT(OUT)     :: qsnow     ! Total H2O snow[kg/m2]
160    REAL, DIMENSION(klon), INTENT(OUT)     :: snowhgt   ! Snow height (m)
161    REAL, DIMENSION(klon), INTENT(OUT)     :: to_ice    ! Snow passed to ice
162    REAL, DIMENSION(klon), INTENT(OUT)     :: sissnow   ! Snow in model (kg/m2)
163                                                                           
164! +--OUTPUT for NetCDF File                                         
165! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                             
166!   REAL, DIMENSION(klonv)                 :: SOsoKL    ! Abs Solar Radiation
167!   REAL, DIMENSION(klonv)                 :: IRsoKL    ! Abs IR    Radiation
168!   REAL, DIMENSION(klonv)                 :: HSsoKL    ! Abs Sensible Ht Flux
169!   REAL, DIMENSION(klonv)                 :: HLsoKL    ! Abs Latent Heat Flux
170!   REAL, DIMENSION(klonv)                 :: HLs_KL    ! Evaporation         
171!   REAL, DIMENSION(klonv)                 :: HLv_KL    ! Transpiration       
172!                             
173!   REAL, DIMENSION(klonv)                 :: SOsoNC    ! Abs Solar Radiation 
174!   REAL, DIMENSION(klonv)                 :: IRsoNC    ! Abs IR    Radiation 
175!   REAL, DIMENSION(klonv)                 :: HSsoNC    ! Abs Sensible Ht Flux
176!   REAL, DIMENSION(klonv)                 :: HLsoNC    ! Abs Latent Heat Flux
177!   REAL, DIMENSION(klonv)                 :: HLs_NC    ! Evaporation     
178!   REAL, DIMENSION(klonv)                 :: HLv_NC    ! Transpiration     
179!                                     
180!   REAL, DIMENSION(klonv,nsoilmx)           :: eta_NC    ! nsoilmx=nsol+1
181!   REAL, DIMENSION(klonv,nsoilmx)           :: tsolNC    !   
182!   REAL, DIMENSION(klonv)                   :: snowNC    !
183!   INTEGER, DIMENSION(klonv)                :: isnoNC    !
184!   INTEGER, DIMENSION(klonv)                :: ispiNC    !     
185!   INTEGER, DIMENSION(klonv)                :: iiceNC    !           
186!   REAL, DIMENSION(klonv)                   :: swaNC     !
187!   REAL, DIMENSION(klonv)                   :: swsNC     !   
188!   INTEGER, DIMENSION(klonv,nsno)           :: istoNC    !
189!   REAL, DIMENSION(klonv,nsno)              :: dzsnNC    !
190!   REAL, DIMENSION(klonv,nsno)              :: rhosnNC   !     
191!   REAL, DIMENSION(klonv,nsno)              :: etasnNC   ! 
192!   REAL, DIMENSION(klonv,nsno)              :: tsnNC     !
193!   REAL, DIMENSION(klonv,nsno)              :: g1snNC    !     
194!   REAL, DIMENSION(klonv,nsno)              :: g2snNC    !                   
195!   REAL, DIMENSION(klonv,nsno)              :: agsnNC    !
196!   REAL, DIMENSION(klonv)                   :: meltNC    !     
197!   REAL, DIMENSION(klonv)                   :: refrNC    !               
198!   REAL, DIMENSION(klonv)                   :: alb1NC    !     
199!   REAL, DIMENSION(klonv)                   :: alb2NC    !         
200!   REAL, DIMENSION(klonv)                   :: rnofNC    !
201                                                                         
202! + Optional Variables:                                                 
203! +--V,  dT(a-s)    Time Moving Averages                                       
204! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~           
205
206! #AW real             V__mem(klon,ntaver)   ! ntaver defined in LMDZ_SL.inc 
207! #AW real             VVmmem(klon)          !                                 
208! #AW common/SVeSBLmem/V__mem,VVmmem         !                                 
209! #AH real             T__mem(klon,ntaver)   !                                 
210! #AH real             dTmmem(klon)          !                                 
211! #AH common/STeSBLmem/T__mem,dTmmem         !                     
212! +--u*, u*T*, u*s* Time Moving Averages                                   
213! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                                       
214! #AM real             u__mem(klon,ntaver)   ! ntaver defined in LMDZ_SL.inc 
215! #AT real             uT_mem(klon,ntaver)   !                                 
216! #AS real             us_mem(klon,ntaver)   !                                 
217! #AM common/S_eSBLmem/u__mem                !                                 
218! #AT.                ,uT_mem                !                                 
219! #AS.                ,us_mem                !
220!! #AH    REAL, DIMENSION(klonv,ntaver)            :: T__mem    !
221!! #AH    REAL, DIMENSION(klonv)                   :: dTmmem    !
222!! #AW    REAL, DIMENSION(klonv,ntaver)            :: V__mem    !     
223!! #AM    REAL, DIMENSION(klonv)                   :: u__mem    !             
224!! #AT    REAL, DIMENSION(klonv)                   :: uT_mem    !
225!! #AS    REAL, DIMENSION(klonv)                   :: us_mem    !     
226!! #AW    REAL, DIMENSION(klonv)                   :: VVmmem    !   
227
228                                                                   
229! +--Internal  Variables                                                       
230! +  ===================                                         
231
232    CHARACTER(len=20)               :: fichnom, fn_outfor ! Name for output file
233    INTEGER                         :: i, ig, ikl, isl, isn, nt
234    INTEGER                         :: gp_outfor, un_outfor
235
236    REAL, PARAMETER                 :: f1=0.5
237    REAL, PARAMETER                 :: sn_upp=5000.,sn_low=500.
238    REAL, PARAMETER                 :: sn_add=400.,sn_div=2.
239                                             ! snow mass upper,lower limit,
240                                             ! added mass/division lowest layer
241    REAL, PARAMETER                 :: c1_zuo=12.960e+4, c2_zuo=2.160e+6
242    REAL, PARAMETER                 :: c3_zuo=1.400e+2,  czemin=1.e-3 
243                                             ! Parameters for drainage
244! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
245! +...        Run Off Parameters                                             
246! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)   
247! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)             
248
249    REAL, DIMENSION(klon)           :: eps0SL          ! surface Emissivity
250    REAL                            :: zsigma, Ua_min, Us_min
251    REAL                            :: lambda          ! Par. soil discret.
252    REAL, DIMENSION(nsoilmx), SAVE  :: dz1,dz2         ! Soil layer thicknesses
253!$OMP THREADPRIVATE(dz1)
254    LOGICAL, SAVE                   :: firstcall=.TRUE.,SnoMod, ok_outfor=.FALSE.!$OMP THREADPRIVATE(firstcall)               
255
256
257
258! + Optional: 
259!c #BW INTEGER                      ::  noUNIT                 
260!c #BW REAL                         ::  BlowST,SnowSB   
261                                         
262                                       
263! +--Internal Variables
264! +  ==================
265
266    INTEGER                         ::  ivg,iso
267
268!========================================================================
269 
270      SnoMod=.true.
271      zsigma=1000.
272      dt__SV=dtime
273      IF (ok_outfor) THEN
274        un_outfor=51                 ! unit number for point output
275        gp_outfor=79! 633 !79        ! grid point number for point output
276        fn_outfor='outfor_SV.dat'
277      END IF
278
279!     write(*,*)'Start of simulation? ',debut        !hj
280      IF (debut) THEN
281        firstcall=.TRUE.
282        INI_SV=.false.
283      ELSE
284        firstcall=.false.
285        INI_SV=.true.
286      END IF
287
288! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
289
290
291! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
292! + INITIALISATION: BEGIN +++
293! + -------------------------
294! +
295! + Compute soil discretization (as for LMDZ)
296! + ----------------------------------------- 
297      IF (firstcall) THEN
298
299! +--Array size
300        klonv=klon
301        knonv=knon
302        write(*,*)'klon',klon,'klonv',klonv,'knon',knon
303
304
305        CALL INIT_VARtSV
306        CALL INIT_VARxSV
307        CALL INIT_VARySV
308
309        eps0SL(:)=0.
310! +--Soil layer thickness                                                   
311! +  ----------------------- 
312!        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
313        CALL get_soil_levels(dz1,dz2,lambda)
314        lambSV=lambda
315        dz1_SV(1:knon,1:) = 0.     
316        dz2_SV(1:knon,1:) = 0.
317                     
318        DO isl =   -nsol,0   
319          dz_dSV(isl) = 0.5e-3*dz2(1-isl)           ! Soil layer thickness
320          DO ikl=1,knon
321            dz1_SV(ikl,isl) = dz1(1-isl)    !1.e-3*
322            dz2_SV(ikl,isl) = dz2(1-isl)    !1.e-3*
323          END DO
324!          IF (knon > 0) THEN
325!            write(*,*)'level:',dz_dSV(isl),dz1_SV(1,isl),dz2_SV(1,isl)
326!          END IF
327        END DO
328
329        DO ikl=1,knon                               
330          eps0SL(ikl )= 1.
331          alb0SV(ikl) = alb_soil(ikl)              ! Soil Albedo       
332
333! + Soil Upward IR Flux, Water Fluxes, roughness length 
334          IRs_SV(ikl) =                               &
335              -eps0SL(ikl)* rsigma*temp_air(ikl)      &  ! Upward IR Flux
336              *temp_air(ikl)*temp_air(ikl)*temp_air(ikl)         
337          TvegSV(ikl) = temp_air(ikl)
338 ! + Soil
339        DO isl =   -nsol,0   
340          TsisSV(ikl,isl) = temp_air(ikl)  !tsoil(ikl,1-isl)  Soil Temperature
341          eta_SV(ikl,isl) = 0.0001         !etasoil(ikl,1-isl)Soil Water[m3/m3]
342          ro__SV(ikl,isl) = 50.      !rosoil(ikl,1-isl) soil water volumic mass
343        END DO           
344        END DO                                           
345
346! +--Surface Fall Line Slope                                                   
347! +  ----------------------- 
348        IF (SnoMod)  THEN               
349          DO ikl=1,knon 
350            slopSV(ikl) = slope(ikl)
351            SWf_SV(ikl) =             &   ! Normalized Decay of the 
352              exp(-dt__SV             &   ! Surficial Water Content 
353              /(c1_zuo                &   !(Zuo and Oerlemans 1996, 
354            +c2_zuo*exp(-c3_zuo*abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
355          END DO                                     
356        END IF           
357                             
358! +--SBL  Characteristics history, not stored at restart.                     
359        DO ikl=1,knon 
360         DO nt=1,ntaver                                       ! #AW
361           V__mem(ikl,nt)=wind_velo(ikl)                      ! #AH
362           T__mem(ikl,nt)=temp_air(ikl)-tsurf(ikl)
363         END DO                                         
364        END DO                                     
365
366
367! + SISVAT_ini (as for use with MAR, but not computing soil layers)
368! + -------------------------------------------------------------
369!        write(*,'(/a)') 'Start SISVAT initialization: SISVAT_ini'
370        CALL SISVAT_ini(knon)
371
372! open output file
373        IF (ok_outfor) THEN
374          open(unit=un_outfor,status='new',file=fn_outfor)         
375          rewind un_outfor     
376          ikl=gp_outfor     ! index sur la grille land ice
377          ig=611            ! index sur la grille globale
378          write(un_outfor,501) fn_outfor, ikl, rlon(ig),rlat(ig) 
379501     format(/,a18,/,'Grid point ',i4,' Long',f9.4,' Lat ',f9.4       &
380     &            ,/,'++++++++++++++++++++++++++++++++++++++++++++++',  &
381     & '++++++++++++++++++++++++++++++++++++++++++++++',                &
382     & '++++++++++++++++++++++++++++++++++++++++++++++',                &
383     &             /,' SWdown + IRdown + Wind   + Temp.  + Humid.   ',  &
384     & '+ Press  +Precip_l+Precip_s+ Tsrf   + Clouds +'                 &
385     & '+ Zenith + BLhgt  + Densair+ Exner +'                           &
386     &            ,/,' sol_SV + IRd_SV + VV__SV + TaT_SV + QaT_SV   ',  &
387     & '+ ps__SV + drr_SV + dsn_SV + Tsf_SV + cld_SV +'                 &
388     & '+ coszSV + za__SV + rhT_SV + ExnrSV+'                           &
389     &            ,/,' W/m2   + W/m2   + m/s    + K      + kg/kg    ',  &
390     & '+ Pa     + kg/m2/s+ kg/m2/s+ K      + /1     +'        &
391     & '+ -      + m      + kg/m3  +       +'                           &
392     &            ,/,'++++++++++++++++++++++++++++++++++++++++++++++',  &
393     & '++++++++++++++++++++++++++++++++++++++++++++++',                &
394     & '++++++++++++++++++++++++++++++++++++++++++++++')
395        END IF
396
397! +--Read restart file
398! +  =================================================   
399       
400! Martin
401       PRINT*, 'On debranche sisvatetat0'
402! Martin
403
404       ! CALL sisvatetat0("startsis.nc",ikl2i)
405
406 
407      END IF  ! firstcall                       
408! +                               
409! +  +++  INITIALISATION:  END  +++                               
410! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
411
412
413
414! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
415! + READ FORCINGS 
416! + ------------------------
417! + Update Forcings for SISVAT given by the LMDZ model.
418! +
419      DO ikl=1,knon
420
421! +--Atmospheric Forcing                                    (INPUT)           
422! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^
423        zSBLSV      = 1000.                         ! [m]               
424        za__SV(ikl) = bl_height(ikl)                ! Height boundary layer [m]
425        Ua_min      = epsi                          !
426! #VM   Ua_min      = 0.2 * sqrt(za__SV(ikl)   )    !                   
427        Ua_min      = 0.2 * sqrt(za__SV(ikl)   )    !                   
428        VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
429        Us_min      = 0.01
430        us__SV(ikl) = max(Us_min, us__SV(ikl) )
431        TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
432        ExnrSV(ikl) = pexner(ikl)                   ! Exner potential         
433        rhT_SV(ikl) = dens_air(ikl)                 ! Air density           
434        QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity
435! #VX   dQa_SV(ikl) = 0. !hj   dtDiff/zsigma(mz)    ! Water Vapor Flux
436        ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
437        p1l_SV(ikl) = p1lay(ikl)                  ! lowest atm. layer press[Pa]
438
439! +--Energy Fluxes                                          (INPUT)           
440! +  ^^^^^^^^^^^^^                                           ^^^^^             
441        coszSV(ikl) = max(czemin,rmu0(ikl))         ! cos(zenith.Dist.) 
442        sol_SV(ikl) = swdown(ikl)                   ! downward Solar 
443        IRd_SV(ikl) = lwdown(ikl)                   ! downward IR   
444        rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.   
445!hj 110511
446!hj 121011        IRs_SV(ikl) = min(IRs_SV(ikl),-epsi)        ! check upward IR   
447!hj
448! +--Water  Fluxes                                          (INPUT)           
449! +  ^^^^^^^^^^^^^                                           ^^^^^             
450        drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
451        dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]
452!c #BS  dbsnow      = -SLussl(i,j,n)                ! Erosion   
453!c #BS.               *dtPhys     *rhT_SV(ikl) /ro_Wat                   
454!c #BS  dsnbSV(ikl) = snow_adv(ikl)  ! min(max(zero,dbsnow)             
455!c #BS.                    /    max(epsi,d_snow),unun)                   
456!c #BS  dbs_SV(ikl) = snow_cont_air(ikl)
457!c #BS                  blowSN(i,j,n)               !          [kg/m2] 
458                                                                             
459! +--Soil/Canopy                                            (INPUT)           
460! +  ^^^^^^^^^^^                                             ^^^^^           
461        alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo
462        AcoHSV(ikl) = AcoefH(ikl) 
463        BcoHSV(ikl) = BcoefH(ikl)                     
464        AcoQSV(ikl) = AcoefQ(ikl) 
465        BcoQSV(ikl) = BcoefQ(ikl)
466        cdH_SV(ikl) = cdragh(ikl)                         
467        Tsf_SV(ikl) = tsurf(ikl)        !hj 12 03 2010
468
469! +--Energy Fluxes                                          (INPUT/OUTPUT)   
470! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^   
471        IF (.not.firstcall) THEN 
472! active hj 110411
473        cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness         
474        END IF
475 
476! +--Time Averages of wind and surface-atmosphere temperature difference
477! +--for turbulence calculations                                       ! #AA
478! Update stored arrays
479      DO nt=1,ntaver-1                                                 ! #AW   
480        V__mem(ikl,nt    ) = V__mem(ikl,nt+1)                          ! #AH   
481        T__mem(ikl,nt    ) = T__mem(ikl,nt+1)                          ! #AA 
482      ENDDO                                                            ! #AW
483        V__mem(ikl,ntaver)=wind_velo(ikl)                              ! #AH
484        T__mem(ikl,ntaver)=temp_air(ikl)-tsurf(ikl)                    ! #AW
485! Calculate averages
486        VVmmem(ikl)        = 0.0                                       ! #AH
487        dTmmem(ikl)        = 0.0                                       ! #AA
488      DO nt=1,ntaver                                                   ! #AW
489        VVmmem(ikl)        = VVmmem(ikl)+V__mem(ikl,nt)                ! #AH
490        dTmmem(ikl)        = dTmmem(ikl)+T__mem(ikl,nt)                ! #AA
491      ENDDO                                                            ! #AW
492        VVmmem(ikl)        = VVmmem(ikl)/ntaver                        ! #AH
493        dTmmem(ikl)        = dTmmem(ikl)/ntaver                   
494      END DO
495
496!                           
497! +  +++  READ FORCINGS:  END  +++   
498! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
499
500! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
501! +  OUTPUT FORCING
502!
503      IF (ok_outfor) THEN
504        ikl=gp_outfor
505        write(un_outfor,5000) sol_SV(ikl), IRd_SV(ikl), VV__SV(ikl),    &   
506     &                        TaT_SV(ikl), QaT_SV(ikl), ps__SV(ikl),    &
507     &                        drr_SV(ikl), dsn_SV(ikl), Tsf_SV(ikl),    &
508     &                        cld_SV(ikl), coszSV(ikl), za__SV(ikl),    &
509     &                        rhT_SV(ikl), ExnrSV(ikl)
510     
511 5000     format(f8.3,' ',f8.3,' ',f8.4,' ',f8.4,' ',f10.8,' ',f8.1,' ', &
512     &           f8.6,' ',f8.6,' ',f8.4,' ',f8.6,' ',f8.6,' ',f8.3,' ',  &
513     &           f8.4,' ',f8.6)
514
515      ENDIF
516!                           
517! +  OUTPUT FORCINGS:  END  +++   
518! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
519
520
521
522! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
523! +--SISVAT EXECUTION                                                         
524! +  ----------------                                                         
525! +             
526!      write(*,*) '  Start SISVAT execution!'
527
528      call  SISVAT(SnoMod,.false.,1)
529!                          BloMod,jjtime)
530
531! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
532! + RETURN RESULTS 
533! + --------------
534! + Return (compressed) SISVAT variables to LMDZ             
535! +
536      DO  ikl=1,knon                  ! use only 1:knon (actual ice sheet..)
537        runoff_lic(ikl)    = RnofSV(ikl)*dtime   ! RunOFF: intensity* time step
538        dflux_s(ikl)       = dSdTSV(ikl)         ! Sens.H.Flux T-Der.
539        dflux_l(ikl)       = dLdTSV(ikl)         ! Latn.H.Flux T-Der.
540        fluxsens(ikl)      = HSs_sv(ikl)         ! HS                 
541        fluxlat(ikl)       = HLs_sv(ikl)         ! HL                 
542        evap(ikl)          = HLs_sv(ikl)/RLVTT   ! Evaporation RLVTT=LhvH2O   
543        z0_new(ikl)        = Z0h_SV(ikl)         ! Moment.Roughn.L.
544
545
546        snow(ikl)          = 0.
547        snowhgt(ikl)       = 0.
548        qsnow(ikl)         = 0.
549        qsol(ikl)          = 0.
550! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
551! +   Check snow thickness, add if too thin, substract if too thick     
552
553        sissnow(ikl)       = 0.  !()
554      DO  isn = 1,isnoSV(ikl)                                               
555        sissnow(ikl)       = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)   
556      END DO
557
558      IF (sissnow(ikl) .LE. sn_low) THEN  !add snow
559      IF (isnoSV(ikl).GE.1) THEN
560        dzsnSV(ikl,1)      = dzsnSV(ikl,1) + sn_add/max(ro__SV(ikl,1),epsi) 
561        toicSV(ikl)        = toicSV(ikl)   - sn_add
562      ELSE
563        write(*,*) 'Attention, bare ice... point ',ikl
564        isnoSV(ikl)        = 1
565        istoSV(ikl,1)      = 100
566!ym        istoSV(ikl)     = 100
567        ro__SV(ikl,1)      = 400.     
568        dzsnSV(ikl,1)      = sn_add/max(ro__SV(ikl,1),epsi)  ! 1.
569        eta_SV(ikl,1)      = epsi
570        TsisSV(ikl,1)      = min(TsisSV(ikl,0),TfSnow-0.2)   
571        G1snSV(ikl,1)      = 99.
572        G2snSV(ikl,1)      = 0.3
573        agsnSV(ikl,1)      = 10.   
574        toicSV(ikl)        = toicSV(ikl)   - sn_add
575      END IF
576      END IF
577
578      IF (sissnow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
579        dzsnSV(ikl,1)      = dzsnSV(ikl,1)/sn_div
580        toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div
581      END IF
582
583        sissnow(ikl)       = 0.  !()
584      DO  isn = 1,isnoSV(ikl)                                               
585        sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)           
586! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
587        snowhgt(ikl) = snowhgt(ikl)+dzsnSV(ikl,isn)             
588        qsnow(ikl)   = qsnow(ikl)+1e03*eta_SV(ikl,isn)*dzsnSV(ikl,isn)   
589      END DO
590        snow(ikl)    = sissnow(ikl)+toicSV(ikl)
591        to_ice(ikl)  = toicSV(ikl)
592! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
593
594
595      DO  isl =   -nsol,0   
596        tsoil(ikl,1-isl)   = TsisSV(ikl,isl)       ! Soil Temperature 
597        qsol(ikl)          = qsol(ikl)                      &   
598                       +eta_SV(ikl,isl) * dz_dSV(isl) 
599      END DO                                               
600        agesno(ikl)        = agsnSV(ikl,isnoSV(ikl))        !          [day]
601
602        alb1(ikl)          = alb1sv(ikl)             ! Albedo VIS 
603        alb2(ikl)          = ((So1dSV-f1)*alb1sv(ikl)                   &
604     &                       +So2dSV*alb2sv(ikl)+So3dSV*alb3sv(ikl))/f1   
605                                                     ! Albedo NIR
606        alb3(ikl)          = alb3sv(ikl)             ! Albedo FIR 
607        tsurf_new(ikl)     = TsfnSV(ikl)
608!hj220711          Tsrfsv(ikl)              ! Surf.Temperature   
609!                TsisSV(ikl,0) *(1-min(1,isnoSV(ikl))) &                   
610!                +TsisSV(ikl,max(1,isnoSV(ikl))) * min(1,isnoSV(ikl))
611!        tsurf_new(ikl)   = max(Ts_Min,tsurf_new(ikl) )
612!        tsurf_new(ikl)   = min(Ts_Max,tsurf_new(ikl) )
613        qsurf(ikl)         = QaT_SV(ikl)
614        emis_new(ikl)      = eps0SL(ikl) 
615
616!!!hjp 230611 sorties
617!!        qsnow(ikl)=TsisSV(ikl,isnoSV(ikl))
618!!        sissnow(ikl)=TsisSV(ikl,isnoSV(ikl)-1)
619!!        snowhgt(ikl)=TsisSV(ikl,isnoSV(ikl)-2)
620!!        qsol(ikl)=dzsnSV(ikl,isnoSV(ikl)-1)
621!!        agesno(ikl)=dzsnSV(ikl,isnoSV(ikl))
622
623
624      END DO
625
626
627! +  -----------------------------                             
628! +  END --- RETURN RESULTS   
629! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
630
631
632! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
633! + Uncompress SISVAT output variables & store
634! +  -----------------------------           
635! +
636!      DO  ikl = 1,knon                                                     
637!        i   = ikl2i(ikl)                             ! Compression index 
638!      IF (i.GE.0) THEN
639!        SOsoNC(i)       = SOsoKL(ikl)                ! Absorb.Sol.Rad.   
640!        IRsoNC(i)       = IRsoKL(ikl)                ! Absorb.IR  Rad.   
641!        HSsoNC(i)       = HSsoKL(ikl)                ! HS                 
642!        HLsoNC(i)       = HLsoKL(ikl)                ! HL                 
643!        HLs_NC(i)       = HLs_KL(ikl)                ! Evaporation         
644!        HLv_NC(i)       = HLv_KL(ikl)                ! Transpiration     
645!       
646!        DO isl =   -nsol,0     
647!          eta_NC(i,1-isl) = eta_SV(ikl,isl)          ! Soil Humidity     
648!          tsolNC(i,1-isl) = TsisSV(ikl,isl)          ! Soil Temperature   
649!        END DO                                                 
650!        snowNC(i)       = snow(ikl)                  ! Snow mass   
651!        isnoNC(i)       = isnoSV(ikl)                ! Nb Snow/Ice Lay.   
652!        ispiNC(i)       = ispiSV(ikl)                ! Nb Supr.Ice Lay.   
653!        iiceNC(i)       = iiceSV(ikl)                ! Nb      Ice Lay.   
654!        swaNC(i)        = rusnSV(ikl)                ! Surficial Water   
655!        swsNC(i)        = SWS_SV(ikl)                ! Surficial Wat.St.   
656!        DO  isn = 1,nsno                                                   
657!          istoNC(i,isn)   = istoSV(ikl,isn)          !            [-]     
658!          dzsnNC(i,isn)   = dzsnSV(ikl,isn)          !            [m]     
659!          rhosnNC(i,isn)  = ro__SV(ikl,isn)          !        [kg/m3]     
660!          etasnNC(i,isn)  = eta_SV(ikl,isn)          !        [m3/m3]     
661!          tsnNC(i,isn)    = TsisSV(ikl,isn)          !            [K]     
662!          g1snNC(i,isn)   = G1snSV(ikl,isn)          ! [-]        [-]     
663!          g2snNC(i,isn)   = G2snSV(ikl,isn)          ! [-] [0.0001 m]     
664!          agsnNC(i,isn)   = agsnSV(ikl,isn)          !          [day]   
665!        END DO                                                           
666!!?c #IB  depsubNC(i)    = wes_SV(ikl)                ! Depo. / Subli.     
667!        meltNC(i)       = wem_SV(ikl)                ! Melting           
668!        refrNC(i)       = wer_SV(ikl)                ! Refreezing         
669!        alb1NC(i)       = alb1sv(ikl)                ! Albedo SW 
670!        alb2NC(i)       = (alb1sv(ikl)+3*alb2sv(ikl)+alb3sv(ikl))/5 
671                                                      ! Albedo LW
672!        rnofNC(i)       = RnofSV(ikl)              ! Run OFF Intensity
673!   
674!!?        SL_z0(i)       = Z0m_SV(ikl)                 ! Moment.Roughn.L.   
675!!?        SL_r0(i,j,n)   = Z0h_SV(ikl)                 ! Heat   Roughn.L. 
676!!!        zWE_NC(i)       = zWE_SV(ikl)                ! Current  *Thick.   
677!!!        zWEcNC(i)       = zWEcSV(ikl)                ! Non-Erod.*Thick.   
678!!!        hSalNC(i)       = hSalSV(ikl)                ! Salt.Layer Height 
679!!!        hsenSL(i) = -SLuts(i,j)  *  rhAir    *cp     ! Sensible Heat Flux
680!!!        hlatSL(i) = -SLuqs(i,j)  *  rhAir    *Lv_H2O ! Latent   Heat Flux
681!      END IF
682!      END DO   
683
684!hj + Put storage to Output file here!                                     
685      IF (lafin) THEN
686
687        fichnom = "restartsis.nc"
688        CALL sisvatredem("restartsis.nc",ikl2i,rlon,rlat)           
689        IF (ok_outfor) THEN
690          close(unit=un_outfor)                   
691        END IF
692      END IF                                                         
693! +  -----------------------------                             
694! +  END --- RETURN RESULTS   
695! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
696
697  END SUBROUTINE surf_sisvat
698
699
700
701
702
703
704
705  SUBROUTINE get_soil_levels(dz1, dz2, lambda)
706! ======================================================================
707! Routine to compute the vertical discretization of the soil in analogy
708! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
709! of SISVAT, therefore it's needed here.
710!
711    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
712    USE mod_phys_lmdz_para
713
714    INCLUDE "dimsoil.h"
715
716    REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
717    REAL, INTENT(OUT)                     :: lambda
718
719
720!-----------------------------------------------------------------------
721!   Depthts:
722!   --------
723    REAL fz,rk,fz1,rk1,rk2
724    REAL min_period, dalph_soil
725    INTEGER ierr,jk
726
727    fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
728
729!    write(*,*)'Start soil level computation'
730!-----------------------------------------------------------------------
731! Calculation of some constants
732! NB! These constants do not depend on the sub-surfaces
733!-----------------------------------------------------------------------
734!-----------------------------------------------------------------------
735!   ground levels
736!   grnd=z/l where l is the skin depth of the diurnal cycle:
737!-----------------------------------------------------------------------
738
739     min_period=1800. ! en secondes
740     dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
741! !$OMP MASTER
742!     IF (is_mpi_root) THEN
743!        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
744!        IF (ierr == 0) THEN ! Read file only if it exists
745!           READ(99,*) min_period
746!           READ(99,*) dalph_soil
747!           PRINT*,'Discretization for the soil model'
748!           PRINT*,'First level e-folding depth',min_period, &
749!                '   dalph',dalph_soil
750!           CLOSE(99)
751!        END IF
752!     ENDIF
753! !$OMP END MASTER
754!     CALL bcast(min_period)
755!     CALL bcast(dalph_soil)
756
757!   la premiere couche represente un dixieme de cycle diurne
758     fz1=SQRT(min_period/3.14)
759     
760     DO jk=1,nsoilmx
761        rk1=jk
762        rk2=jk-1
763        dz2(jk)=fz(rk1)-fz(rk2)
764     ENDDO
765     DO jk=1,nsoilmx-1
766        rk1=jk+.5
767        rk2=jk-.5
768        dz1(jk)=1./(fz(rk1)-fz(rk2))
769     ENDDO
770     lambda=fz(.5)*dz1(1)
771     PRINT*,'full layers, intermediate layers (seconds)'
772     DO jk=1,nsoilmx
773        rk=jk
774        rk1=jk+.5
775        rk2=jk-.5
776        PRINT *,'fz=', &
777             fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
778     ENDDO
779
780  END SUBROUTINE get_soil_levels
781
782  SUBROUTINE SISVAT_ini(knon)                                                     
783                                                                               
784!C +------------------------------------------------------------------------+ 
785!C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR |
786!C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters | 
787!C +------------------------------------------------------------------------+ 
788!C |   PARAMETERS:  klonv: Total Number of columns =                        | 
789!C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       | 
790!C |                     X       Number of Mosaic Cell per grid box         | 
791!C |                                                                        | 
792!C |   INPUT:   dt__SV   : Time  Step                                   [s] |
793!C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] | 
794!C |                                                                        | 
795!C |   OUTPUT:  RF__SV   : Root Fraction in Layer isl                   [-] | 
796!C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] | 
797!C |            etamSV   : Soil Minimum Humidity                    [m3/m3] | 
798!C |                      (based on a prescribed Soil Relative Humidity)    | 
799!C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       | 
800!C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        | 
801!C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   | 
802!C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] | 
803!C |            dzsnSV(0): Soil first Layer Thickness                   [m] | 
804!C |            dzmiSV   : Distance between two contiguous levels       [m] | 
805!C |            dz78SV   : 7/8 (Layer Thickness)                        [m] |
806!C |            dz34SV   : 3/4 (Layer Thickness)                        [m] |
807!C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] |
808!C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] |
809!C |            dtz_SV   : dt/dz                                      [s/m] |
810!C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
811!C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      | 
812!C |            Explic   : Explicit Parameter = 1.0 - Implic                | 
813!C |                                                                        |
814!C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
815!C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
816!C | #          #SH: Hapex-Sahel Values                                     ! 
817!C |                                                                        |
818!C +------------------------------------------------------------------------+ 
819!                                                                             
820!                                                                             
821                                                                             
822!C +--Global Variables                                                         
823!C +  ================         
824
825      USE VARphy, ra_earth=>ra                                           
826      USE VAR_SV                                                     
827      USE VARdSV                                                         
828      USE VAR0SV                                                           
829      USE VARxSV
830      USE VARtSV
831      USE VARxSV
832      USE VARySV
833      IMPLICIT NONE                                                           
834                                                                             
835                                                                             
836                                                                             
837!C +--Arguments                                                     
838!C +  ==================                                                     
839       INTEGER,INTENT(IN) ::  knon                                       
840
841!C +--Internal Variables                                                     
842!C +  ==================                                                     
843                                                                               
844      INTEGER ::  ivt   ,ist   ,ivg   ,ikl   ,isl   ,isn   ,ikh               
845      INTEGER ::  misl_2,nisl_2                                             
846      REAL    ::  zDepth                                                       
847      REAL    ::  d__eta,eta__1,eta__2,Khyd_1,Khyd_2                           
848      REAL,PARAMETER  ::  RHsMin=  0.001        ! Min.Soil Relative Humidity   
849      REAL    ::  PsiMax                        ! Max.Soil Water    Potential
850      REAL    ::  a_Khyd,b_Khyd                 ! Piecewis.Water Conductivity
851
852
853!c #WR REAL    ::  Khyd_x,Khyd_y                                               
854                                                                             
855                                                                             
856                                                                 
857!C +--Non Time Dependant SISVAT parameters                                   
858!C +  ====================================                               
859                                                                             
860!C +--Soil Discretization                                                     
861!C +  -------------------                                                     
862                                                                             
863!C +--Numerical Scheme Parameters                                             
864!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^                                             
865        Implic = 0.75                           ! 0.5  <==> Crank-Nicholson 
866        Explic = 1.00 - Implic                  !                           
867                                                                             
868!C +--Soil/Snow Layers Indices                                               
869!C +  ^^^^^^^^^^^^^^^^^^^^^^^^                                               
870      DO  isl=-nsol,0                                                         
871        islpSV(isl) =           isl+1                                         
872        islpSV(isl) = min(      islpSV(isl),0)                               
873        islmSV(isl) =           isl-1                                         
874        islmSV(isl) = max(-nsol,islmSV(isl))                                 
875      END DO                                                                 
876                                                                               
877      DO  isn=1,nsno                                                           
878        isnpSV(isn) =           isn+1                                         
879        isnpSV(isn) = min(      isnpSV(isn),nsno)                           
880      END DO                                                                 
881                                                                             
882!C +--Soil      Layers Thicknesses                                             
883!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
884! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!   
885!c #kd IF (nsol.gt.4)                                              THEN       
886!c #kd   DO isl=-5,-nsol,-1                                                   
887!c #kd     dz_dSV(isl)=   1.                                                 
888!c #kd   END DO                                                               
889!c #kd END IF                                                                 
890!                                                                             
891!      IF (nsol.ne.4)                                              THEN       
892!        DO isl= 0,-nsol,-1                                                   
893!          misl_2 =     -mod(isl,2)                                         
894!          nisl_2 =         -isl/2                                           
895!          dz_dSV(isl)=(((1-misl_2) * 0.001                                   
896!     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.             
897!C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm           
898!                                                                             
899!c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001                                   
900!c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.             
901!                                                                             
902!c #05     dz_dSV(isl)=(((1-misl_2) * 0.001                                   
903!c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5             
904!        END DO                                                               
905!          dz_dSV(0)  =               0.001                                   
906!          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004         
907!      END IF                                                                 
908                                                                               
909        zz_dSV      = 0.                                                     
910      DO  isl=-nsol,0                                                       
911        dzmiSV(isl) = 0.500*(dz_dSV(isl)        +dz_dSV(islmSV(isl)))       
912        dziiSV(isl) = 0.500* dz_dSV(isl)        /dzmiSV(isl)                 
913        dzi_SV(isl) = 0.500* dz_dSV(islmSV(isl))/dzmiSV(isl)                   
914        dtz_SV(isl) =        dt__SV             /dz_dSV(isl)                 
915        dz78SV(isl) = 0.875* dz_dSV(isl)                                     
916        dz34SV(isl) = 0.750* dz_dSV(isl)                                     
917        dz_8SV(isl) = 0.125* dz_dSV(isl)                                       
918        dzAvSV(isl) = 0.125* dz_dSV(islmSV(isl))                        &
919     &              + 0.750* dz_dSV(isl)                                &
920     &              + 0.125* dz_dSV(islpSV(isl))                         
921!c #ER   dz78SV(isl) =        dz_dSV(isl)                                     
922!c #ER   dz34SV(isl) =        dz_dSV(isl)                                     
923!c #ER   dz_8SV(isl) = 0.                                                     
924!c #ER   dzAvSV(isl) =        dz_dSV(isl)                                     
925        zz_dSV      = zz_dSV+dz_dSV(isl)                                     
926      END DO                                                                 
927      DO ikl=1,knon !v                                                         
928        dzsnSV(ikl,0) =      dz_dSV(0)                                       
929      END DO                                                                 
930                                                                             
931!C +--Conversion to a 50 m Swab Ocean Discretization                           
932!C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                           
933        OcndSV = 0.                                                           
934      DO isl=-nsol,0                                                         
935        OcndSV = OcndSV +dz_dSV(isl)                                           
936      END DO                                                                   
937        OcndSV = 50.    /OcndSV                                               
938                                                                               
939                                                                               
940!C +--Secondary Vegetation Parameters                                         
941!C +  -------------------------------                                         
942!                                                                             
943!C +--Minimum Stomatal Resistance (Hapex Sahel Data)                           
944!C +  (Taylor et al. 1997, J.Hydrol 188-189, p.1047)                           
945!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                           
946!c #SH DO ivg=1,3                       !                                     
947!c #SH   StodSV(ivg) = 210.             ! Millet                               
948!c #SH END DO                           !                                     
949!c #SH   StodSV(  4) = 120.             ! Sparse Tiger Bush                   
950!c #SH DO ivg=5,6                       !                                     
951!c #SH   StodSV(ivg) =  80.             ! Dense  Tiger Bush                   
952!c #SH END DO                           !                                     
953!c #SH   StodSV(  7) =  80.             ! Low    Trees (Fallow)               
954!c #SH   StodSV( 10) =  80.             !                                     
955!                                                                             
956!C +--Minimum Stomatal Resistance (Tropical Forest)                         
957!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                           
958!c #SH   StodSV(  8) =  60.             ! Medium Trees                       
959!c #SH   StodSV( 11) =  60.             !                                     
960!c #SH   StodSV(  9) =  40.             ! High   Trees                       
961!c #SH   StodSV( 12) =  40.             !                                     
962                                                                             
963!C +--Root Fraction                                                           
964!C +  ^^^^^^^^^^^^^                                                         
965!C +  * GENERAL REFERENCE                                                     
966!C +    Jackson et al., 1996: A global analysis of root distributions for     
967!C +    terrestrial biomes. In Oecologia, 108, 389-411.                       
968                                                                             
969!C +  * ROOT PROFILE                                                           
970!C +    The cumulative root fraction Y is given by                             
971!C +        Y = 1 - beta**d   with d    the depth (in cm),                   
972!C +                               beta a coefficient (vegetation dependent).
973                                                                             
974!C +  * BETA VALUES (for 11 world biomes)                                     
975!C +  1  boreal forest                0.943                                   
976!C +  2  crops                        0.961                                   
977!C +  3  desert                       0.975                                   
978!C +  4  sclerophyllous shrubs        0.964                                   
979!C +  5  temperate coniferous forest  0.976                                   
980!C +  6  temperate deciduous forest   0.966                                   
981!C +  7  temperate grassland          0.943                                   
982!C +  8  tropical deciduous forest    0.961                                   
983!C +  9  tropical evergreen forest    0.962                                   
984!C +  10 tropical grassland savanna   0.972                                   
985!C +  11 tundra                       0.914                                   
986!                                                                             
987!C +  * ADVISED BETA VALUES FOR MAR                                           
988!C +    (see 'block data SISVAT_dat', variable rbtdSV)                         
989!C +                                                                           
990!C +    SVAT veg. type         default      West Africa                       
991!C +    0  barren soil         0.000        0.000                             
992!C +    1  crops low           0.961 (2)    0.961 (2)                         
993!C +    2  crops medium        0.961 (2)    0.961 (2)                         
994!C +    3  crops high          0.961 (2)    0.961 (2)                         
995!C +    4  grass low           0.943 (7)    0.943 (7)                         
996!C +    5  grass medium        0.943 (7)    0.964 (4)                         
997!C +    6  grass high          0.943 (7)    0.972 (10)                       
998!C +    7  broadleaf low       0.966 (6)    0.968 (4,10)                     
999!C +    8  broadleaf medium    0.966 (6)    0.962 (8,9)                       
1000!C +    9  broadleaf high      0.966 (6)    0.962 (8,9)                       
1001!C +    10 needleleaf low      0.976 (5)    0.971 (5,6)                       
1002!C +    11 needleleaf medium   0.976 (5)    0.976 (5)                         
1003!C +    12 needleleaf high     0.976 (5)    0.976 (5)                         
1004                                                                             
1005!C +    Numbers between brackets refer to Jackson's biomes. For more details 
1006!C +    about some choices, see the correspondance between the IGBP and SVAT   
1007!C +    vegetation classes (i.e. in NESTOR).                                 
1008                                                                             
1009!C +  * WARNING                                                               
1010!C +    Most of the roots are located in the first 2 m of soil. The root       
1011!C +    fraction per layer depends on the definition of the soil layer         
1012!C +    thickness. It will get wrong if a thick layer is defined around 2 m   
1013!C +    deep.                                                                 
1014                                                                               
1015      write(*,'(/a)') 'ROOT PROFILES (Jackson, 1996) :'                       
1016                                                                             
1017      DO ivt = 0, nvgt                                                       
1018        zDepth = 0.                                                           
1019        DO isl = 0, -nsol, -1                                                 
1020          IF (ivt .ne. 0) THEN                                               
1021            RF__SV(ivt,isl) =  rbtdSV(ivt)**zDepth *                     &
1022     &                         (1. - rbtdSV(ivt)**(dz_dSV(isl)*100) )       
1023            zDepth = zDepth + dz_dSV(isl)*100  !in cm                       
1024          ELSE                                                                 
1025            RF__SV(ivt,isl) = 0.                                             
1026          END IF                                                             
1027        END DO                                                               
1028        write(*,'(a,i2,a,i3,a,99f10.5:)')                                &
1029     &       '  RF__SV(', ivt, ',', -nsol, ':0) =', RF__SV(ivt,:)           
1030      END DO                                                                 
1031!      write(6, format(                                                  &
1032!     &  '  NOTE: If root fraction is not close to 0  around 2 m deep,', &
1033!     &/,'        Then you should redefine the soil layer thicknesses.', &
1034!     &/,'        See the code for more details.'))                         
1035                                                                             
1036                                                                               
1037!C +--Secondary Soil       Parameters                                         
1038!C +  -------------------------------                                         
1039                                                                               
1040      DO  ist=0,nsot                                                           
1041         rocsSV(ist)=(1.0-etadSV(ist))*1.2E+6  ! Soil Contrib. to (ro c)_s   
1042         s1__SV(ist)=     bCHdSV(ist)          & ! Factor of (eta)**(b+2)     
1043     &  *psidSV(ist)     *Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)     
1044     & /(etadSV(ist)**(   bCHdSV(ist)+3.))     !                             
1045         s2__SV(ist)=     Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)     
1046     & /(etadSV(ist)**(2.*bCHdSV(ist)+3.))     !    in DR97, Eqn.(3.35)       
1047                                                                             
1048!C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)     
1049!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^   
1050         Psimax = -(log(RHsMin))/7.2E-5        ! DR97, Eqn 3.15 Inversion   
1051         etamSV(ist) =  etadSV(ist)                                      &
1052     &         *(PsiMax/psidSV(ist))**(-min(10.,1./bCHdSV(ist)))             
1053      END DO                                                                 
1054         etamSV(12)  =  0.                                                   
1055                                                                             
1056!C +--Piecewise Hydraulic Conductivity Profiles                               
1057!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
1058      DO   ist=0,nsot                                                         
1059!c #WR      write(6, format(' Type |    etaSat | No |    eta__1 |    eta__2 |', &
1060!c #WR&           '    Khyd_1 |    Khyd_x |    Khyd_2 |    Khyd_y |'        &
1061!c #WR&         /,' -----+-----------+----+-----------+-----------+',       &
1062!c #WR&           '-----------+-----------+-----------+-----------+'))       
1063                                                                             
1064          d__eta          =  etadSV(ist)/nkhy                                 
1065          eta__1          =  0.                                             
1066          eta__2          =  d__eta                                           
1067        DO ikh=0,nkhy                                                         
1068          Khyd_1          =  s2__SV(ist)             & ! DR97, Eqn.(3.35)     
1069     &  *(eta__1      **(2. *bCHdSV(ist)+3.))        !                       
1070          Khyd_2          =  s2__SV(ist)             &!                       
1071     &  *(eta__2      **(2. *bCHdSV(ist)+3.))        !                       
1072                                                                             
1073          a_Khyd          = (Khyd_2-Khyd_1)/d__eta   !                       
1074          b_Khyd          =  Khyd_1-a_Khyd *eta__1   !                       
1075!c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !                       
1076!c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !                       
1077          aKdtSV(ist,ikh) =  a_Khyd       * dt__SV   !                       
1078          bKdtSV(ist,ikh) =  b_Khyd       * dt__SV   !                         
1079!c #WR     write(6,format(i5,' |',e10.2,' |',i3,' |', 6(e10.2,' |'))      &
1080!c #WR&           ist,etadSV(ist),ikh,eta__1,                             &
1081!c #WR&          eta__2,Khyd_1,Khyd_x,Khyd_2,Khyd_y                       
1082                                                                         
1083          eta__1          = eta__1  + d__eta                             
1084          eta__2          = eta__2  + d__eta                             
1085        END DO                                                             
1086      END DO                                                               
1087
1088
1089!c
1090! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1091! + INITIALISATION: ARRAYS on grid  +++
1092! + -----------------------------------
1093! +                                                     
1094        DO ikl=1,knon                               
1095
1096! + Water Fluxes, roughness length 
1097          us__SV(ikl) = 0.25                             ! Frict. Velocity 
1098                                      ! (approx ETH-Camp Lefebre etal CD 2005
1099          uts_SV(ikl) = 1.                               ! u*T*  arbitrary 
1100          uqs_SV(ikl) = 0.25 ! turb_vel(ikl,3)           ! u*q*    "
1101          uss_SV(ikl) = 0.25 ! turb_vel(ikl,4)           ! u*s*    "
1102!!c #AE   usthSV(ikl) = 0.
1103          Z0h_SV(ikl) = 0. ! rlength(ikl,9)
1104          Z0m_SV(ikl) = 0.                               ! Moment.Roughn.L. 
1105          Z0h_SV(ikl) = 0.                               ! Heat   Roughn.L. 
1106!!c #OR   Z0roSV(ikl) = 0.                               ! Orogr. Roughn.L.
1107          LMO_SV(ikl) = 0.
1108          VVs_SV(ikl) = 0.
1109          RRs_SV(ikl) = 0.
1110          DDs_SV(ikl) = 0.
1111
1112! + Snow
1113          isnoSV(ikl) = 0                        ! # snow layers
1114          BufsSV(ikl) = 0.   
1115          zzsnsv(ikl,:) = 0.                         ! Snow pack thickness
1116          qsnoSV(ikl) = 0.                        ! BL snow content 
1117          zWEcSV(ikl) = 0.
1118          dbs_SV(ikl) = 0.
1119          dsnbSV(ikl) = 0.
1120          esnbSV(ikl) = 0.
1121          BrosSV(ikl) = 0.
1122          BG1sSV(ikl) = 0.         
1123          BG2sSV(ikl) = 0.
1124          SWS_SV(ikl) = 0.
1125          RnofSV(ikl) = 0.                             ! RunOFF Intensity
1126          toicSV(ikl) = 0.         !hj preli
1127! + Clouds !hj1602
1128          cld_SV(ikl) = 0.           
1129
1130! + Vegetation
1131          LSmask(ikl) = 1                          ! Land/Sea   Mask   
1132          isotSV(ikl) = 12                         ! Soil       Type   
1133          iWaFSV(ikl) = 1                          ! Soil Drainage         
1134          Rootsv(ikl,:) = 0.
1135          ivgtSV(ikl) = 0                          ! Vegetation Type   
1136          LAI0SV(ikl) = 0                          ! LAI               
1137          glf0SV(ikl) = 0                          ! Green Leaf Frac.
1138!          TVegSV(ikl) = 0.
1139          SnCaSV(ikl) = 0.
1140          rrCaSV(ikl) = 0.
1141          pSivSV(ikl) = 0.
1142
1143! + z0(Sastrugi)                                       
1144! #SZ     Z0SaSL(ikl) = 0.                                           
1145! #ZM   DO nt=1,ntavSL                                                     
1146! #ZM     SLn_z0(ikl,nt) = 0.5e-6                                       
1147! #ZM     SLn_b0(ikl,nt) = 0.5e-6                                       
1148! #ZM     SLn_r0(ikl,nt) = 0.5e-6                                       
1149! #ZM   END DO
1150
1151
1152! +--z0(Orography Roughness)                                 
1153! #OR     SL_z0 (ikl) = min(SL_z0 (ikl),zsigma/3.)                   
1154! #OR     SLzoro(ikl) = min(SLzoro(ikl),zsigma/3.)                   
1155
1156! +--SBL  Characteristics                                                         
1157! #AA   DO nt=  1,ntaver                                                     
1158! #AW     V_0aSL(ikl,nt)  = ssvSL(ikl)   
1159! #AH     dT0aSL(ikl,nt)  = temp_air(ikl)-tsurf(ikl)                     
1160! #AA   END DO                                                               
1161
1162
1163        END DO
1164
1165
1166                                                                           
1167      return                                                               
1168
1169  END SUBROUTINE SISVAT_ini
1170
1171
1172
1173
1174
1175
1176
1177!***************************************************************************
1178
1179    SUBROUTINE sisvatetat0 (fichnom,ikl2i)
1180
1181    USE dimphy
1182    USE mod_grid_phy_lmdz
1183    USE mod_phys_lmdz_para
1184
1185    USE iostart
1186    USE VAR_SV
1187    USE VARdSV
1188    USE VARxSV         
1189    USE VARtSV
1190    USE indice_sol_mod
1191
1192      IMPLICIT none
1193!======================================================================
1194! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1195! Objet: Ecriture de l'etat de redemarrage pour SISVAT
1196!======================================================================
1197    include "netcdf.inc"
1198!    include "indicesol.h"
1199
1200    include "dimsoil.h"
1201    include "clesphys.h"
1202    include "thermcell.h"
1203    include "compbl.h"
1204
1205!======================================================================
1206    CHARACTER(LEN=*) :: fichnom
1207
1208
1209    INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1210    REAL, DIMENSION(klon) :: rlon
1211    REAL, DIMENSION(klon) :: rlat
1212
1213! les variables globales ecrites dans le fichier restart
1214    REAL, DIMENSION(klon) :: isno
1215    REAL, DIMENSION(klon) :: ispi
1216    REAL, DIMENSION(klon) :: iice
1217    REAL, DIMENSION(klon) :: rusn
1218    REAL, DIMENSION(klon, nsno) :: isto
1219
1220    REAL, DIMENSION(klon, nsismx) :: Tsis
1221    REAL, DIMENSION(klon, nsismx) :: eta
1222    REAL, DIMENSION(klon, nsismx) :: ro
1223
1224    REAL, DIMENSION(klon, nsno) :: dzsn     
1225    REAL, DIMENSION(klon, nsno) :: G1sn
1226    REAL, DIMENSION(klon, nsno) :: G2sn
1227    REAL, DIMENSION(klon, nsno) :: agsn
1228
1229    REAL, DIMENSION(klon) :: toic
1230
1231!    REAL, DIMENSION(klon) :: IRs
1232!    REAL, DIMENSION(klon) :: LMO
1233!    REAL, DIMENSION(klon) :: Bufs
1234!    REAL, DIMENSION(klon, 9) :: rlength
1235!    REAL, DIMENSION(klon, 5) :: turb_vel
1236
1237    INTEGER  :: isl, ikl, i, isn , errT, erreta, errro, errdz, snopts
1238    CHARACTER (len=2) :: str2
1239    LOGICAL :: found
1240 
1241    errT=0
1242    errro=0
1243    erreta=0
1244    errdz=0
1245    snopts=0
1246! Ouvrir le fichier contenant l'etat initial:
1247
1248    CALL open_startphy(fichnom)
1249
1250! Lecture des latitudes, longitudes (coordonnees):
1251
1252      CALL get_field("latitude",rlat,found)
1253      CALL get_field("longitude",rlon,found)
1254
1255      CALL get_field("n_snows", isno,found)
1256      IF (.NOT. found) THEN
1257        PRINT*, 'phyetat0: Le champ <n_snows> est absent'
1258        PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
1259      ENDIF
1260
1261      CALL get_field("n_ice_top",ispi,found)
1262      CALL get_field("n_ice",iice,found)
1263      CALL get_field("surf_water",rusn,found)
1264!      IF (.NOT. found) THEN
1265!        PRINT*, 'phyetat0: Le champ <surf_water> est absent'
1266!        rusn(:)=0. 
1267!      ENDIF
1268
1269
1270      CALL get_field("to_ice",toic,found)
1271      IF (.NOT. found) THEN
1272        PRINT*, 'phyetat0: Le champ <to_ice> est absent'
1273        toic(:)=0. 
1274      ENDIF
1275
1276
1277!      CALL get_field("IR_soil",IRs,found)
1278!      CALL get_field("LMO",LMO,found)
1279!      CALL get_field("snow_buffer",Bufs,found)
1280!      DO i = 1, 5
1281!            WRITE(str2,'(i2.2)') i
1282!            CALL get_field("turb_veloc"//str2, &
1283!                          turb_vel(:,i),found)
1284!      ENDDO
1285!      DO i = 1, 9
1286!            WRITE(str2,'(i2.2)') i
1287!            CALL get_field("rough_length"//str2, &
1288!                          rlength(:,i),found)
1289!      ENDDO
1290      DO isn = 1,nsno
1291        IF (isn.LE.99) THEN
1292            WRITE(str2,'(i2.2)') isn
1293            CALL get_field("AGESNOW"//str2, &
1294                          agsn(:,isn),found)
1295        ELSE
1296            PRINT*, "Trop de couches"
1297            CALL abort
1298        ENDIF
1299      ENDDO
1300      DO isn = 1,nsno
1301        IF (isn.LE.99) THEN
1302            WRITE(str2,'(i2.2)') isn
1303            CALL get_field("DZSNOW"//str2, &
1304                          dzsn(:,isn),found)
1305        ELSE
1306            PRINT*, "Trop de couches"
1307            CALL abort
1308        ENDIF
1309      ENDDO
1310      DO isn = 1,nsno
1311        IF (isn.LE.99) THEN
1312            WRITE(str2,'(i2.2)') isn
1313            CALL get_field("G2SNOW"//str2, &
1314                          G2sn(:,isn),found)
1315        ELSE
1316            PRINT*, "Trop de couches"
1317            CALL abort
1318        ENDIF
1319      ENDDO
1320      DO isn = 1,nsno
1321        IF (isn.LE.99) THEN
1322            WRITE(str2,'(i2.2)') isn
1323            CALL get_field("G1SNOW"//str2, &
1324                          G1sn(:,isn),found)
1325        ELSE
1326            PRINT*, "Trop de couches"
1327            CALL abort
1328        ENDIF
1329      ENDDO
1330      DO isn = 1,nsismx
1331        IF (isn.LE.99) THEN
1332            WRITE(str2,'(i2.2)') isn
1333            CALL get_field("ETA"//str2, &
1334                          eta(:,isn),found)
1335        ELSE
1336            PRINT*, "Trop de couches"
1337            CALL abort
1338        ENDIF
1339      ENDDO
1340      DO isn = 1,nsismx
1341        IF (isn.LE.99) THEN
1342            WRITE(str2,'(i2.2)') isn
1343            CALL get_field("RO"//str2, &
1344                          ro(:,isn),found)
1345        ELSE
1346            PRINT*, "Trop de couches"
1347            CALL abort
1348        ENDIF
1349      ENDDO
1350      DO isn = 1,nsismx
1351        IF (isn.LE.99) THEN
1352            WRITE(str2,'(i2.2)') isn
1353            CALL get_field("TSS"//str2, &
1354                          Tsis(:,isn),found)
1355        ELSE
1356            PRINT*, "Trop de couches"
1357            CALL abort
1358        ENDIF
1359      ENDDO
1360      DO isn = 1,nsno
1361        IF (isn.LE.99) THEN
1362            WRITE(str2,'(i2.2)') isn
1363            CALL get_field("HISTORY"//str2, &
1364                          isto(:,isn),found)
1365        ELSE
1366            PRINT*, "Trop de couches"
1367            CALL abort
1368        ENDIF
1369      ENDDO
1370      write(*,*)'Read ',fichnom,' finished!!'
1371
1372!*********************************************************************************
1373! Compress restart file variables for SISVAT
1374
1375
1376      isnoSV(:)        = 0       
1377      ispiSV(:)        = 0
1378      iiceSV(:)        = 0                               
1379       
1380      eta_SV(:,1:nsno) = 0.     
1381      TsisSV(:,1:nsno) = 0.         
1382      istoSV(:,:)      = 0
1383      ro__SV(:,1:nsno) = 0.       
1384      G1snSV(:,:)      = 0.   
1385      G2snSV(:,:)      = 0.
1386      dzsnSV(:,1:nsno) = 0.       
1387      agsnSV(:,:)      = 0.
1388      rusnSV(:)        = 0.   
1389      toicSV(:)        = 0.
1390
1391      DO  ikl = 1,klon                                                   
1392          i   = ikl2i(ikl)   
1393          IF (i > 0) THEN
1394              isnoSV(ikl)     = INT(isno(i))          ! Nb Snow/Ice Lay.   
1395              ispiSV(ikl)     = INT(ispi(i))          ! Nb Supr.Ice Lay. 
1396              iiceSV(ikl)     = INT(iice(i))          ! Nb      Ice Lay.   
1397                                                                             
1398!!              IRs_SV(ikl)     = IRs(i)
1399!              LMO_SV(ikl)     = LMO(i)               !?             
1400!              us__SV(ikl)     = turb_vel(i,1)       
1401!              uts_SV(ikl)     = turb_vel(i,2)
1402!              uqs_SV(ikl)     = turb_vel(i,3)
1403!              uss_SV(ikl)     = turb_vel(i,4)
1404!              usthSV(ikl)     = turb_vel(i,5)
1405!              Z0m_SV(ikl)     = rlength(i,1)         ! Moment.Roughn.L.
1406!              Z0mmSV(ikl)     = rlength(i,2)         ! Moment.Roughn.L.   
1407!              Z0mnSV(ikl)     = rlength(i,3)         ! Moment.Roughn.L.   
1408!              Z0SaSV(ikl)     = rlength(i,4)         ! Moment.Roughn.L.   
1409!              Z0e_SV(ikl)     = rlength(i,5)         ! Moment.Roughn.L.   
1410!              Z0emSV(ikl)     = rlength(i,6)         ! Moment.Roughn.L.   
1411!              Z0enSV(ikl)     = rlength(i,7)         ! Moment.Roughn.L.   
1412!              Z0roSV(ikl)     = rlength(i,8)         ! Moment.Roughn.L.   
1413!!              Z0h_SV(ikl)     = rlength(i,9)         ! Moment.Roughn.L.   
1414
1415            DO isl =   -nsol,0   
1416              ro__SV(ikl,isl) = ro(i,nsno+1-isl)       !                   
1417              eta_SV(ikl,isl) = eta(i,nsno+1-isl)         ! Soil Humidity     
1418!hjp 15/10/2010
1419              IF (eta_SV(ikl,isl) <= 1.e-6) THEN          !hj check
1420                eta_SV(ikl,isl) = 1.e-6
1421              ENDIF
1422              TsisSV(ikl,isl) = Tsis(i,nsno+1-isl)        ! Soil Temperature 
1423              IF (TsisSV(ikl,isl) <= 1.) THEN          !hj check
1424!                errT=errT+1
1425                TsisSV(ikl,isl) = 273.15
1426              ENDIF
1427
1428            END DO       
1429            write(*,*)'Copy histo', ikl
1430
1431 
1432              istoSV(ikl,0) = 0                     ! Snow     History
1433
1434              G1snSV(ikl,0) = 0.5                   !      [-]     
1435              G2snSV(ikl,0) = 3.                    ! [-] [0.0001 m]
1436              dzsnSV(ikl,0) = dz_dSV(0)             !            [m]
1437              agsnSV(ikl,0) = 0.                    !          [day]     
1438   
1439            DO  isn = 1,isnoSV(ikl) !nsno     
1440              snopts=snopts+1
1441              IF (isto(i,isn) > 10.) THEN          !hj check
1442                write(*,*)'Irregular isto',ikl,i,isn,isto(i,isn)
1443                isto(i,isn) = 1.
1444              ENDIF
1445
1446              istoSV(ikl,isn) = INT(isto(i,isn))     ! Snow     History
1447              ro__SV(ikl,isn) = ro(i,isn)            !        [kg/m3]     
1448              eta_SV(ikl,isn) = eta(i,isn)           !        [m3/m3] 
1449              TsisSV(ikl,isn) = Tsis(i,isn)          !            [K] 
1450
1451             IF (TsisSV(ikl,isn) <= 1.) THEN          !hj check
1452              errT=errT+1
1453              TsisSV(ikl,isn) = TsisSV(ikl,0)
1454             ENDIF 
1455             IF (TsisSV(ikl,isn) <= 1.) THEN          !hj check
1456              TsisSV(ikl,isn) = 263.15
1457             ENDIF
1458             IF (eta_SV(ikl,isn) < 1.e-9) THEN          !hj check
1459              eta_SV(ikl,isn) = 1.e-6 
1460              erreta=erreta+1
1461             ENDIF 
1462             IF (ro__SV(ikl,isn) <= 10.) THEN          !hj check
1463              ro__SV(ikl,isn) = 11.
1464              errro=errro+1
1465             ENDIF
1466              write(*,*)ikl,i,isn,Tsis(i,isn),G1sn(i,isn)
1467              G1snSV(ikl,isn) = G1sn(i,isn)          ! [-]        [-]     
1468              G2snSV(ikl,isn) = G2sn(i,isn)          ! [-] [0.0001 m]
1469              dzsnSV(ikl,isn) = dzsn(i,isn)          !            [m]
1470!             IF (dzsnSV(ikl,isn) < 5.e-6) THEN          !hj check
1471!              dzsnSV(ikl,isn) = 0.000005
1472!              errdz=errdz+1
1473!             ENDIF         
1474              agsnSV(ikl,isn) = agsn(i,isn)          !          [day]     
1475            END DO 
1476              rusnSV(ikl)     = rusn(i)              ! Surficial Water   
1477              toicSV(ikl)     = toic(i)              ! bilan snow to ice   
1478!!              BufsSV(ikl)     = Bufs(i)   
1479          END IF                       
1480        END DO   
1481!        write(*,*)snopts,'snow pts',errT,' T errors,',errro,'ro,'
1482!        write(*,*)'   ',erreta,'eta,',errdz,'dz'
1483    END SUBROUTINE sisvatetat0
1484
1485
1486
1487
1488    SUBROUTINE sisvatredem (fichnom,ikl2i,rlon,rlat)
1489!======================================================================
1490! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1491! Objet: Ecriture de l'etat de redemarrage pour SISVAT
1492!======================================================================
1493    USE mod_grid_phy_lmdz
1494    USE mod_phys_lmdz_para
1495    USE iostart
1496    USE VAR_SV
1497    USE VARxSV         
1498    USE VARySV !hj tmp 12 03 2010
1499    USE VARtSV
1500    USE indice_sol_mod
1501
1502    IMPLICIT none
1503
1504    include "netcdf.inc"
1505!    include "indicesol.h"
1506    include "dimsoil.h"
1507    include "clesphys.h"
1508    include "thermcell.h"
1509    include "compbl.h"
1510
1511!======================================================================
1512
1513    CHARACTER(LEN=*) :: fichnom
1514    INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1515    REAL, DIMENSION(klon), INTENT(IN) :: rlon
1516    REAL, DIMENSION(klon), INTENT(IN) :: rlat
1517
1518! les variables globales ecrites dans le fichier restart
1519    REAL, DIMENSION(klon) :: isno
1520    REAL, DIMENSION(klon) :: ispi
1521    REAL, DIMENSION(klon) :: iice
1522    REAL, DIMENSION(klon, nsnowmx) :: isto
1523
1524    REAL, DIMENSION(klon, nsismx) :: Tsis
1525    REAL, DIMENSION(klon, nsismx) :: eta
1526    REAL, DIMENSION(klon, nsnowmx) :: dzsn
1527    REAL, DIMENSION(klon, nsismx) :: ro       
1528    REAL, DIMENSION(klon, nsnowmx) :: G1sn
1529    REAL, DIMENSION(klon, nsnowmx) :: G2sn
1530    REAL, DIMENSION(klon, nsnowmx) :: agsn
1531    REAL, DIMENSION(klon) :: IRs
1532    REAL, DIMENSION(klon) :: LMO
1533    REAL, DIMENSION(klon) :: rusn
1534    REAL, DIMENSION(klon) :: toic
1535    REAL, DIMENSION(klon) :: Bufs
1536    REAL, DIMENSION(klon) :: alb1,alb2,alb3
1537    REAL, DIMENSION(klon, 9) :: rlength
1538    REAL, DIMENSION(klon, 5) :: turb_vel
1539
1540    INTEGER isl, ikl, i, isn
1541    CHARACTER (len=2) :: str2
[3518]1542    INTEGER           :: pass
[1865]1543
1544      isno(:)       = 0       
1545      ispi(:)       = 0
1546      iice(:)       = 0                               
1547      IRs(:)        = 0.
1548      LMO(:)        = 0.                             
1549      turb_vel(:,:) = 0.
1550      rlength(:,:)  = 0.                                 
1551      eta(:,:)      = 0.     
1552      Tsis(:,:)     = 0.         
1553      isto(:,:)     = 0
1554      ro(:,:)       = 0.       
1555      G1sn(:,:)     = 0.   
1556      G2sn(:,:)     = 0.
1557      dzsn(:,:)     = 0.       
1558      agsn(:,:)     = 0.
1559      rusn(:)       = 0.   
1560      toic(:)       = 0.   
1561      Bufs(:)       = 0.   
1562      alb1(:)       = 0.
1563      alb2(:)       = 0.
1564      alb3(:)       = 0.
1565
1566!***************************************************************************
1567! Uncompress SISVAT output variables for storage
1568           
1569      DO  ikl = 1,klon                                                   
1570        i   = ikl2i(ikl)
1571      IF (i > 0) THEN
1572        isno(i)       = 1.*isnoSV(ikl)               ! Nb Snow/Ice Lay.   
1573        ispi(i)       = 1.*ispiSV(ikl)               ! Nb Supr.Ice Lay.   
1574        iice(i)       = 1.*iiceSV(ikl)               ! Nb      Ice Lay.       
1575   
1576!        IRs(i)        = IRs_SV(ikl)
1577!        LMO(i)        = LMO_SV(ikl)                             
1578!        turb_vel(i,1) = us__SV(ikl)
1579!        turb_vel(i,2) = uts_SV(ikl)
1580!        turb_vel(i,3) = uqs_SV(ikl)
1581!        turb_vel(i,4) = uss_SV(ikl)
1582!        turb_vel(i,5) = usthSV(ikl)
1583!        rlength(i,1)  = Z0m_SV(ikl)                  ! Moment.Roughn.L.
1584!        rlength(i,2)  = Z0mmSV(ikl)                  ! 
1585!        rlength(i,3)  = Z0mnSV(ikl)                  ! 
1586!        rlength(i,4)  = Z0SaSV(ikl)                  !   
1587!        rlength(i,5)  = Z0e_SV(ikl)                  !   
1588!        rlength(i,6)  = Z0emSV(ikl)                  !   
1589!        rlength(i,7)  = Z0enSV(ikl)                  !   
1590!        rlength(i,8)  = Z0roSV(ikl)                  !   
1591!        rlength(i,9)  = Z0h_SV(ikl)                  !   
1592
1593        DO isl =   -nsol,0                           !                   
1594          eta(i,nsno+1-isl)  = eta_SV(ikl,isl)            ! Soil Humidity     
1595          Tsis(i,nsno+1-isl) = TsisSV(ikl,isl)            ! Soil Temperature   
1596          ro(i,nsno+1-isl)   = ro__SV(ikl,isl)            !        [kg/m3]   
1597        END DO       
1598 
1599 
1600        DO  isn = 1,nsno             
1601          isto(i,isn)   = 1.*istoSV(ikl,isn)         ! Snow     History
1602          ro(i,isn)     = ro__SV(ikl,isn)            !        [kg/m3]     
1603          eta(i,isn)    = eta_SV(ikl,isn)            !        [m3/m3] 
1604          Tsis(i,isn)   = TsisSV(ikl,isn)            !            [K]   
1605          G1sn(i,isn)   = G1snSV(ikl,isn)            ! [-]        [-]     
1606          G2sn(i,isn)   = G2snSV(ikl,isn)            ! [-] [0.0001 m]
1607          dzsn(i,isn)   = dzsnSV(ikl,isn)            !            [m]         
1608          agsn(i,isn)   = agsnSV(ikl,isn)            !          [day]     
1609        END DO 
1610        rusn(i)       = rusnSV(ikl)                  ! Surficial Water 
1611        toic(i)       = toicSV(ikl)                  ! Surficial Water 
1612        alb1(i)       = alb1sv(ikl)
1613        alb2(i)       = alb2sv(ikl)
1614        alb3(i)       = alb3sv(ikl)
1615!        Bufs(i)       = BufsSV(ikl)     
1616      END IF                     
1617      END DO                                               
1618
1619      CALL open_restartphy(fichnom)
[3518]1620      DO pass = 1, 2
1621        CALL put_field(pass,"longitude", &
[1865]1622                    "Longitudes de la grille physique",rlon)     
[3518]1623        CALL put_field(pass,"latitude","Latitudes de la grille physique",rlat)
1624 
1625        CALL put_field(pass,"n_snows", "number of snow/ice layers",isno)
1626        CALL put_field(pass,"n_ice_top", "number of top ice layers",ispi)
1627        CALL put_field(pass,"n_ice", "number of ice layers",iice)
1628        CALL put_field(pass,"IR_soil", "Soil IR flux",IRs)
1629        CALL put_field(pass,"LMO", "Monin-Obukhov Scale",LMO)
1630        CALL put_field(pass,"surf_water", "Surficial water",rusn)
1631        CALL put_field(pass,"snow_buffer", "Snow buffer layer",Bufs)
1632        CALL put_field(pass,"alb_1", "albedo sw",alb1)
1633        CALL put_field(pass,"alb_2", "albedo nIR",alb2)
1634        CALL put_field(pass,"alb_3", "albedo fIR",alb3)
1635        CALL put_field(pass,"to_ice", "Snow passed to ice",toic)
[1865]1636
1637 !     DO i = 1, 5
1638 !       WRITE(str2,'(i2.2)') i
[3518]1639 !       CALL put_field(pass,"turb_veloc"//str2, &
[1865]1640 !                      "various turbulent velocities"//str2, &
1641 !                      turb_vel(:,i))
1642 !     ENDDO
1643 !     DO i = 1, 9
1644 !       WRITE(str2,'(i2.2)') i
[3518]1645 !       CALL put_field(pass,"rough_length"//str2, &
[1865]1646 !                      "various roughness lengths"//str2, &
1647 !                      rlength(:,i))
1648 !     ENDDO
[3518]1649        DO isn = 1,nsno
1650          IF (isn.LE.99) THEN
1651            WRITE(str2,'(i2.2)') isn
1652            CALL put_field(pass,"AGESNOW"//str2, &
[1865]1653                         "Age de la neige layer No."//str2, &
1654                         agsn(:,isn))
[3518]1655          ELSE
1656            PRINT*, "Trop de couches"
1657            CALL abort
1658          ENDIF
1659        ENDDO
1660        DO isn = 1,nsno
1661          IF (isn.LE.99) THEN
1662            WRITE(str2,'(i2.2)') isn
1663            CALL put_field(pass,"DZSNOW"//str2, &
[1865]1664                         "Snow/ice thickness layer No."//str2, &
1665                         dzsn(:,isn))
[3518]1666          ELSE
1667            PRINT*, "Trop de couches"
1668            CALL abort
1669          ENDIF
1670        ENDDO
1671        DO isn = 1,nsno
1672          IF (isn.LE.99) THEN
1673            WRITE(str2,'(i2.2)') isn
1674            CALL put_field(pass,"G2SNOW"//str2, &
[1865]1675                         "Snow Property 2, layer No."//str2, &
1676                         G2sn(:,isn))
[3518]1677          ELSE
1678            PRINT*, "Trop de couches"
1679            CALL abort
1680          ENDIF
1681        ENDDO
1682        DO isn = 1,nsno
1683          IF (isn.LE.99) THEN
1684            WRITE(str2,'(i2.2)') isn
1685            CALL put_field(pass,"G1SNOW"//str2, &
[1865]1686                         "Snow Property 1, layer No."//str2, &
1687                         G1sn(:,isn))
[3518]1688          ELSE
1689            PRINT*, "Trop de couches"
1690            CALL abort
1691          ENDIF
1692        ENDDO
1693        DO isn = 1,nsismx
1694          IF (isn.LE.99) THEN
1695            WRITE(str2,'(i2.2)') isn
1696            CALL put_field(pass,"ETA"//str2, &
[1865]1697                         "Soil/snow water content layer No."//str2, &
1698                         eta(:,isn))
[3518]1699          ELSE
1700            PRINT*, "Trop de couches"
1701            CALL abort
1702          ENDIF
1703        ENDDO
1704        DO isn = 1,nsismx   !nsno
1705          IF (isn.LE.99) THEN
1706            WRITE(str2,'(i2.2)') isn
1707            CALL put_field(pass,"RO"//str2, &
1708                           "Snow density layer No."//str2, &
1709                           ro(:,isn))
1710          ELSE
1711            PRINT*, "Trop de couches"
1712            CALL abort
1713          ENDIF
1714        ENDDO
1715        DO isn = 1,nsismx
1716          IF (isn.LE.99) THEN
1717            WRITE(str2,'(i2.2)') isn
1718            CALL put_field(pass,"TSS"//str2, &
1719                           "Soil/snow temperature layer No."//str2, &
1720                           Tsis(:,isn))
1721          ELSE
1722            PRINT*, "Trop de couches"
1723            CALL abort
1724          ENDIF
1725        ENDDO
1726        DO isn = 1,nsno
1727          IF (isn.LE.99) THEN
1728            WRITE(str2,'(i2.2)') isn
1729            CALL put_field(pass,"HISTORY"//str2, &
1730                           "Snow history layer No."//str2, &
1731                           isto(:,isn))
1732          ELSE
1733            PRINT*, "Trop de couches"
1734            CALL abort
1735          ENDIF
1736        ENDDO
[1865]1737      ENDDO
1738
1739  END SUBROUTINE sisvatredem
1740
1741END MODULE surf_sisvat_mod
Note: See TracBrowser for help on using the repository browser.