Ignore:
Timestamp:
Nov 19, 2021, 4:58:59 PM (3 years ago)
Author:
lguez
Message:

Sync latest trunk changes to Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
1 deleted
11 edited
1 copied

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/VARphy.F90

    r3792 r4013  
    2626      INTEGER, PARAMETER ::  iun=1                                             
    2727      REAL, PARAMETER    ::  zer0 = 0.0e+0, half = 0.5e+0, un_1 = 1.0e+0,     &
    28      &                       eps6 = 1.0e-6, R_1000=1.e3      
     28     &                       eps6 = 1.0e-6, R_1000=1.e3   
    2929      REAL, PARAMETER    ::  zero = 0.0e+0, demi = 0.5e+0, unun = 1.0e+0,     &
    3030     &                       epsi = 1.0e-6, eps9 = 1.0e-9         
    … …  
    9191! A1.6 Turbulent and molecular diffusion
    9292!----------------------------------------
    93       REAL, PARAMETER    ::  A_MolV = 1.35e-5, vonKrm = 0.40e0
     93      REAL, PARAMETER    ::  A_MolV = 1.35e-5, vonKrm = 0.40e0, r_turb=3.0
     94      REAL, PARAMETER    ::  A_turb=5.8, akmol=1.35e-5
    9495!C +...                A_MolV: Air Viscosity                 = 1.35d-5 m2/s   
    9596!C +                   vonKrm: von Karman constant           = 0.4           
    96                                                                          
     97!C +                   r_turb:   Turbulent Diffusivities Ratio K*/Km       
     98!C +                   A_turb:    Stability  Coefficient Moment                                 
     99!C +                   Air Viscosity                 = 1.35d-5 m2/s                 
     100
     101
    97102
    98103END MODULE VARphy

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/VARtSV.F90

    r3792 r4013  
    4141
    4242  SUBROUTINE INIT_VARtSV
     43 
    4344  IMPLICIT NONE
    44  
     45
     46  INTEGER ikl
     47
     48
     49
     50
     51
     52
     53
    4554      ALLOCATE(toicSV(klonv))
    4655
    … …  
    5968      ALLOCATE(rsolSV(klonv))                ! Radiation balance surface
    6069
     70      DO ikl=1,klonv       
     71           
     72         toicSV(ikl)   = 0.
     73         dz1_SV(ikl,:) = 0.
     74         dz2_SV(ikl,:) = 0.
     75         Tsf_SV(ikl)   = 0.
     76         TsfnSV(ikl)   = 0.
     77         AcoHSV(ikl)   = 0.
     78         BcoHSV(ikl)   = 0.
     79         AcoQSV(ikl)   = 0.
     80         ps__SV(ikl)   = 0.
     81         p1l_SV(ikl)   = 0.
     82         cdH_SV(ikl)   = 0.
     83         cdM_SV(ikl)   = 0.
     84         rsolSV(ikl)   = 0.
     85      END DO
     86
     87
     88
    6189  END SUBROUTINE INIT_VARtSV
    6290

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/VARxSV.F90

    r3792 r4013  
    6767      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QaT_SV  ! SBL Top   Specific Humidity     
    6868!$OMP THREADPRIVATE(QaT_SV)
     69      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QsT_SV  ! SBL Top   Specific Humidity
     70!$OMP THREADPRIVATE(QsT_SV)
    6971      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dQa_SV  ! SBL Flux  Limitation of Qa     
    7072!$OMP THREADPRIVATE(dQa_SV)
    … …  
    7880
    7981                                                             
    80       REAL,SAVE                      ::   zSBLSV  ! SBL Height (Initial Value)     
    81 !$OMP THREADPRIVATE(zSBLSV)
    8282      REAL,SAVE                      ::   dt__SV  ! Time Step                       
    8383!$OMP THREADPRIVATE(dt__SV)
    … …  
    160160      REAL,ALLOCATABLE,SAVE    ::   agsnSV(:,:)  ! Snow Age                       
    161161!$OMP THREADPRIVATE(agsnSV)
     162      REAL,ALLOCATABLE,SAVE    ::   DOPsnSV(:,:)  ! Snow optical diameter [m]
     163!$OMP THREADPRIVATE(DOPsnSV)
    162164      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BufsSV  ! Snow Buffer Layer               
    163165!$OMP THREADPRIVATE(BufsSV)
    … …  
    260262      ALLOCATE(dLdTSV(klonv))  ! Latent   Heat Flux T Derivat.   
    261263      ALLOCATE(rhT_SV(klonv))  ! SBL Top   Air  Density         
    262       ALLOCATE(QaT_SV(klonv))  ! SBL Top   Specific Humidity     
     264      ALLOCATE(QaT_SV(klonv))  ! SBL Top   Specific Humidity   
     265      ALLOCATE(QsT_SV(klonv))  ! surface   Specific Humidity
    263266      ALLOCATE(dQa_SV(klonv))  ! SBL Flux  Limitation of Qa     
    264267      ALLOCATE(qsnoSV(klonv))  ! SBL Mean  Snow       Content   
    … …  
    309312      ALLOCATE(dzsnSV(klonv,    0:nsno))  ! Snow Layer  Thickness           
    310313      ALLOCATE(agsnSV(klonv,    0:nsno))  ! Snow Age                       
     314      ALLOCATE(DOPsnSV(klonv,    0:nsno))  ! Snow Optical diameter                       
    311315      ALLOCATE(BufsSV(klonv))  ! Snow Buffer Layer               
    312316      ALLOCATE(rusnSV(klonv))  ! Surficial   Water               
    … …  
    339343
    340344      DO ikl=1,klonv       
    341         LSmask(ikl)   = 0
    342         isotSV(ikl)   = 0               
    343         iWaFSV(ikl)   = 0
    344         isnoSV(ikl)   = 0 
    345         ispiSV(ikl)   = 0
    346         iiceSV(ikl)   = 0         
    347         istoSV(ikl,:) = 0
    348         ii__SV(ikl)   = 0
    349         jj__SV(ikl)   = 0
    350         nn__SV(ikl)   = 0
     345
     346
     347      isnoSV(ikl)  =0.       
     348      ispiSV(ikl)  =0.   
     349      iiceSV(ikl)  =0.       
     350      istoSV(ikl,:)=0.                                                                       
     351      alb_SV(ikl)  =0.     
     352      emi_SV(ikl)  =0.   
     353      IRs_SV(ikl)  =0.     
     354      LMO_SV(ikl)  =0.
     355      us__SV(ikl)  =0.
     356      uts_SV(ikl)  =0.
     357      cutsSV(ikl)  =0.
     358      uqs_SV(ikl)  =0.
     359      uss_SV(ikl)  =0.
     360      usthSV(ikl)  =0.
     361      rCDmSV(ikl)  =0.
     362      rCDhSV(ikl)  =0.
     363      Z0m_SV(ikl)  =0.
     364      Z0mmSV(ikl)  =0.
     365      Z0mnSV(ikl)  =0.
     366      Z0roSV(ikl)  =0.
     367      Z0SaSV(ikl)  =0.
     368      Z0e_SV(ikl)  =0.
     369      Z0emSV(ikl)  =0.
     370      Z0enSV(ikl)  =0.
     371      Z0h_SV(ikl)  =0.
     372      Z0hmSV(ikl)  =0.
     373      Z0hnSV(ikl)  =0.
     374                                                                         
     375                                                                       
     376      TsisSV(ikl,:)  =0.
     377      ro__SV(ikl,:)  =0.
     378      eta_SV(ikl,:)  =0.
     379      G1snSV(ikl,:)  =0. 
     380      G2snSV(ikl,:)  =0.         
     381      dzsnSV(ikl,:)  =0.     
     382      agsnSV(ikl,:)  =0.             
     383      DOPsnSV(ikl,:) =0.                 
     384      BufsSV(ikl)  =0.         
     385      rusnSV(ikl)  =0.     
     386      SWf_SV(ikl)  =0.       
     387      SWS_SV(ikl)  =0.
     388      HFraSV(ikl)  =0.       
     389                                                                           
     390      zWE_SV(ikl)  =0.
     391      zWEcSV(ikl)  =0.
     392      wem_SV(ikl)  =0.
     393      wer_SV(ikl)  =0.
     394      wes_SV(ikl)  =0.
     395      zn4_SV(ikl)  =0.
     396      zn5_SV(ikl)  =0.                                                       
     397                                       
     398                                                                                                                                                         
     399      ii__SV(ikl)  =0.
     400      jj__SV(ikl)  =0.
     401      nn__SV(ikl)  =0.
     402                                                                               
     403      IRu_SV(ikl)  =0.
     404      hSalSV(ikl)  =0.     
     405      qSalSV(ikl)  =0.
     406      RnofSV(ikl)  =0.   
     407      RuofSV(ikl,:)  =0.
     408
     409
     410
     411
    351412      END DO
    352413  END SUBROUTINE INIT_VARxSV

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/VARySV.F90

    r3792 r4013  
    2222      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb3sv  ! Surface Albedo FIR     
    2323!$OMP THREADPRIVATE(alb3sv)
    24 
     24      REAL, DIMENSION(:,:),SAVE,ALLOCATABLE  ::   alb6sv  ! 6 band-albedo   
     25!$OMP THREADPRIVATE(alb6sv)
    2526      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   albssv  ! Soil               Albedo [-]   
    2627!$OMP THREADPRIVATE(albssv)
    … …  
    8384      ALLOCATE(alb2sv(klonv))  ! Surface Albedo NIR     
    8485      ALLOCATE(alb3sv(klonv))  ! Surface Albedo FIR       
     86      ALLOCATE(alb6sv(klonv,6))! 6-band  Albedo     
    8587
    8688      !
    … …  
    110112
    111113      DO ikl=1,klonv
    112         NLaysv(ikl)   = 0
    113         i_thin(ikl)   = 0
    114         LIndsv(ikl)   = 0
     114
     115      NLaysv(ikl) =0.
     116      i_thin(ikl) =0.
     117      LIndsv(ikl) =0.
     118      albisv(ikl) =0.
     119      alb1sv(ikl) =0.
     120      alb2sv(ikl) =0.   
     121      alb3sv(ikl) =0.
     122      alb6sv(ikl,:)=0.
     123      albssv(ikl) =0.
     124      SoSosv(ikl) =0.
     125      Eso_sv(ikl) =0.
     126      HSv_sv(ikl) =0.   
     127      HLv_sv(ikl) =0.
     128      HSs_sv(ikl) =0.       
     129      HLs_sv(ikl) =0.
     130      sqrCm0(ikl) =0.   
     131      sqrCh0(ikl) =0.   
     132      Lx_H2O(ikl) =0.
     133      ram_sv(ikl) =0.
     134      rah_sv(ikl) =0.
     135      Fh__sv(ikl) =0.     
     136      dFh_sv(ikl) =0.
     137      Evp_sv(ikl) =0.
     138      EvT_sv(ikl) =0.
     139      LSdzsv(ikl) =0.
     140      Tsrfsv(ikl) =0.
     141      sEX_sv(ikl,:)  =0.
     142      zzsnsv(ikl,:)  =0.
     143      psi_sv(ikl,:)  =0.
     144      Khydsv(ikl,:)  =0.
     145      EExcsv(ikl)  =0.   
     146
     147
    115148      END DO
    116149

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/inlandsis.F

    r3792 r4013  
    1       subroutine INLANDSIS(SnoMod,BloMod,jjtime)
     1      subroutine INLANDSIS(SnoMod,BloMod,jjtime,debut)
    22
    33      USE dimphy
    … …  
    173173      USE VARySV
    174174      USE VARtSV
    175       USE surface_data, only: iflag_tsurf_inlandsis
    176 
     175      USE surface_data, ONLY: is_ok_z0h_rn,
     176     .                        is_ok_density_kotlyakov,
     177     .                        prescribed_z0m_snow,
     178     .                        iflag_z0m_snow,
     179     .                        iflag_tsurf_inlandsis,
     180     .                        iflag_temp_inlandsis,
     181     .                        discret_xf, buf_sph_pol,buf_siz_pol     
    177182
    178183      IMPLICIT NONE
    … …  
    180185      logical   SnoMod
    181186      logical   BloMod
     187      logical   debut
    182188      integer   jjtime
    183189
    … …  
    213219      integer   IceMsk,IcIndx(klonv)          !      Ice / No      Ice Mask
    214220      integer   SnoMsk                        ! Snow     / No Snow     Mask
    215 
    216221      real      roSMin,roSMax,roSn_1,roSn_2,roSn_3   ! Fallen Snow Density (PAHAUT)
    217222      real      Dendr1,Dendr2,Dendr3          ! Fallen Snow Dendric.(GIRAUD)
    218223      real      Spher1,Spher2,Spher3,Spher4   ! Fallen Snow Spheric.(GIRAUD)
    219224      real      Polair                        ! Polar  Snow Switch
    220       real      PorSno,Por_BS,Salt_f,PorRef   !
     225      real      PorSno,Salt_f,PorRef   !
    221226c #sw real      PorVol,rWater                 !
    222227c #sw real      rusNEW,rdzNEW,etaNEW          !
    … …  
    244249      real      Z0m_Sn,Z0m_90                 ! Snow  Surface Roughness Length
    245250      real      SnoWat                        ! Snow Layer    Switch
    246 c #RN real      rstar,alors                   !
    247 c #RN real      rstar0,rstar1,rstar2          !
     251      real      rstar,alors                   !
     252      real      rstar0,rstar1,rstar2          !
    248253      real      SameOK                        ! 1. => Same Type of Grains
    249254      real      G1same                        ! Averaged G1,  same Grains
    … …  
    263268      real      Sph_av                        ! Averaged    Grain Spher.
    264269      real      Den_av                        ! Averaged    Grain Dendr.
    265       real      DendOK                        ! 1. => Average is  Dendr.
    266270      real      G1diff                        ! Averaged G1, diff. Grains
    267271      real      G2diff                        ! Averaged G2, diff. Grains
    … …  
    277281      real      tt_c,vv_c                     ! Critical param.
    278282      real      tt_tmp,vv_tmp,vv_virt         ! Temporary variables
    279       logical   density_kotlyakov             ! .true. if Kotlyakov 1961
    280283      real      e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations
    281284      real      zm1, zm2, coefslope                    ! variables for surface temperature extrapolation
    282 
     285! for Aeolian erosion and blowing snow
     286      integer   nit   ,iit
     287      real      Fac                           ! Correc. factor for drift ratio
     288      real      dusuth,signus
     289      real      sss__F,sss__N
     290      real      sss__K,sss__G
     291      real      us_127,us_227,us_327,us_427,us_527
     292      real      VVa_OK, usuth0
     293      real      ssstar
     294      real      SblPom
     295      real      rCd10n                        ! Square root of drag coefficient
     296      real      DendOK                        ! Dendricity Switch
     297      real      SaltOK                        ! Saltation  Switch
     298      real      MeltOK                        ! Saltation  Switch (Melting Snow)
     299      real      SnowOK                        ! Pack Top   Switch
     300      real      SaltM1,SaltM2,SaltMo,SaltMx   ! Saltation  Parameters
     301      real      ShearX, ShearS                ! Arg. Max Shear Stress
     302      real      Por_BS                        ! Snow Porosity
     303      real      Salt_us                       ! New thresh.friction velocity u*t
     304      real      Fac_Mo,ArguSi,FacRho          ! Numerical factors for u*t
     305      real      SaltSI(klonv,0:nsno)          ! Snow Drift Index              !
     306      real      MIN_Mo                        ! Minimum Mobility Fresh Fallen *
     307      character*3    qsalt_param              ! Switch for saltation flux param.
     308      character*3    usth_param               ! Switch for u*t param
    283309
    284310
    … …  
    287313
    288314      data      T__Min / 200.00/              ! Minimum realistic Temperature
    289       data      TaPole / 263.15/              ! Maximum Polar     Temperature
    290       data      roSMin /  30.  /              ! Minimum Snow  Density
     315      data      TaPole / 268.15/              ! Maximum Polar Temperature (value from C. Agosta)
     316      data      roSMin / 300.  /              ! Minimum Snow  Density
    291317      data      roSMax / 400.  /              ! Max Fresh Snow Density
    292318      data      tt_c   / -2.0  /              ! Critical Temp. (degC)
    … …  
    305331      data      EmiWat /   0.99999999/        ! Emissivity of a Water Area
    306332      data      EmiSno /   0.99999999/        ! Emissivity of Snow
     333
    307334     
    308335!     DATA      Emissivities                  ! Pielke, 1984, pp. 383,409
    … …  
    321348      data      Z0_ICE/    0.0010/            ! Sea-Ice Z0 = 0.0010 m (Andreas)
    322349!                                             !    (Ice Station Weddel -- ISW)
     350! for aerolian erosion
     351      data      SblPom/ 1.27/   ! Lower Boundary Height Parameter
     352C +                             ! for Suspension
     353C +                             ! Pommeroy, Gray and Landine 1993,
     354C +                             ! J. Hydrology, 144(8) p.169
     355      data      nit   / 5   /   ! us(is0,uth) recursivity: Nb Iterations
     356cc#AE data      qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p
     357      data      qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray
     358cc#AE data      usth_param/"lis"/  ! u*t from Liston et al. 2007
     359      data      usth_param/"gal"/  ! u*t from Gallee et al. 2001
     360      data      SaltMx/-5.83e-2/
     361
    323362      vk2    =  vonKrm  *  vonKrm             ! Square of Von Karman Constant
    324363
    … …  
    352391
    353392
    354 ! Blowing Particles Threshold Friction velocity
    355 ! =============================================
    356 
    357 c #AE       usthSV(ikl) =                     1.0e+2
    358 !          END DO
    359 !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
    360 
    361 
    362 
    363 
    364 ! Contribution of Snow to the Surface Snow Pack
    365 ! =============================================
    366 
    367       IF (SnoMod)                                                 THEN
    368 
    369 
    370  
    371 C +--Blowing Snow
    372 C +  ------------
    373  
    374         IF (BloMod) then
    375          if (klonv.eq.1) then
     393
     394
     395
     396      IF (SnoMod)                            THEN
     397
     398 
     399C +--Aeolian erosion and Blowing Snow
     400C +==================================
     401
     402
     403
     404        DO ikl=1,knonv
     405            usthSV(ikl) =                     1.0e+2
     406        END DO
     407
     408
     409        IF (BloMod) THEN
     410 
     411        if (klonv.eq.1) then
    376412          if(isnoSV(1).ge.2                   .and.
    377      .       TsisSV(1,max(1,isnoSV(1)))<273.  .and.
    378      .       ro__SV(1,max(1,isnoSV(1)))<500.  .and.
    379      .       eta_SV(1,max(1,isnoSV(1)))<epsi) then
     413     .         TsisSV(1,max(1,isnoSV(1)))<273.  .and.
     414     .         ro__SV(1,max(1,isnoSV(1)))<500.  .and.
     415     .         eta_SV(1,max(1,isnoSV(1)))<epsi) then
    380416C +                       **********
    381417                     call SISVAT_BSn
    … …  
    384420                     call SISVAT_BSn
    385421C +                       **********
    386          endif
    387         ENDIF
    388  
    389 
    390 
     422        endif
     423
     424
     425
     426
     427
     428
     429
     430! Calculate threshold erosion velocity for next time step
     431! Unlike in sisvat, computation is of threshold velocity made here (instead of sisvaesbl)
     432! since we do not use sisvatesbl for the coupling with LMDZ
     433
     434C +--Computation of threshold friction velocity for snow erosion
     435C ---------------------------------------------------------------
     436
     437        rCd10n =  1. / 26.5 ! Vt / u*t = 26.5
     438                     ! Budd et al. 1965, Antarct. Res. Series Fig.13
     439                     ! ratio developped during assumed neutral conditions
     440 
     441
     442C +--Snow Properties
     443C +  ~~~~~~~~~~~~~~~
     444
     445        DO ikl = 1,knonv
     446
     447          isn      =  isnoSV(ikl)
     448
     449
     450 
     451          DendOK   =  max(zero,sign(unun,epsi-G1snSV(ikl,isn)  ))  !
     452          SaltOK   =  min(1   , max(istdSV(2)-istoSV(ikl,isn),0))  !
     453          MeltOK   =     (unun                                     !
     454     .             -max(zero,sign(unun,TfSnow-epsi                 !
     455     .             -TsisSV(ikl,isn)  )))                           ! Melting Snow
     456     .             *  min(unun,DendOK                              !
     457     .                  +(1.-DendOK)                               !
     458     .                      *sign(unun,     G2snSV(ikl,isn)-1.0))  ! 1.0 for 1mm
     459          SnowOK   =  min(1   , max(isnoSV(ikl)      +1 -isn ,0))  ! Snow Switch
     460 
     461          G1snSV(ikl,isn) =      SnowOK *    G1snSV(ikl,isn)
     462     .                  + (1.- SnowOK)*min(G1snSV(ikl,isn),G1_dSV)
     463          G2snSV(ikl,isn) =      SnowOK *    G2snSV(ikl,isn)
     464     .                  + (1.- SnowOK)*min(G2snSV(ikl,isn),G1_dSV)
     465 
     466          SaltOK   =  min(unun, SaltOK + MeltOK) * SnowOK
     467 
     468 
     469C +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
     470C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     471          SaltM1   = -0.750e-2 * G1snSV(ikl,isn)
     472     .             -0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case
     473C +     CAUTION:  Guyomarc'h & Merindol Dendricity Sign is +
     474C +     ^^^^^^^^                    MAR Dendricity Sign is -
     475          SaltM2   = -0.833d-2 * G1snSV(ikl,isn)
     476     .             -0.583d-2 * G2snSV(ikl,isn)+ 0.833d00 !non-dendritic case
     477 
     478c       SaltMo   = (DendOK   * SaltM1 + (1.-DendOK) *     SaltM2       )
     479          SaltMo   = 0.625 !SaltMo pour d=s=0.5
     480 
     481!weighting SaltMo with surface snow density (Vionnet et al. 2012)
     482cc#AE   FacRho   = 1.25 - 0.0042 * ro__SV(ikl,isn)
     483cc#AE   SaltMo   = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow
     484          MIN_Mo   =  0.
     485c       SaltMo   =  max(SaltMo,MIN_Mo)
     486c       SaltMo   =  SaltOK   * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx)
     487c #TUNE SaltMo   =  SaltOK   * SaltMo - (1.-SaltOK) *     0.9500
     488          SaltMo   =  max(SaltMo,epsi-unun)
     489 
     490C +--Influence of Density on Threshold Shear Stress
     491C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     492          Por_BS =  1. - 300. / ro_Ice
     493          ShearS = Por_BS / (1.-Por_BS)
     494C +...         SheaBS =  Arg(sqrt(shear = max shear stress in snow)):
     495C +            shear  =  3.420d00 * exp(-(Por_BS      +Por_BS)
     496C +  .                                  /(unun        -Por_BS))
     497C +            SheaBS :  see de Montmollin         (1978),
     498C +                      These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124
     499 
     500C +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
     501C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     502          ArguSi      =     -0.085 *us__SV(ikl)/rCd10n
     503!V=u*/sqrt(CD) eqs 2 to 4 Gallee et al. 2001
     504 
     505          SaltSI(ikl,isn) = -2.868 * exp(ArguSi) + 1 + SaltMo
     506 
     507
     508C +--Threshold Friction Velocity
     509C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     510          if(ro__SV(ikl,isn)>300.) then
     511             Por_BS      =  1.000       - ro__SV(ikl,isn)     /ro_Ice
     512          else
     513             Por_BS      =  1.000  - 300. /ro_Ice
     514          endif
     515 
     516          ShearX =  Por_BS/max(epsi,1.-Por_BS)
     517          Fac_Mo = exp(-ShearX+ShearS)
     518C +     Gallee et al., 2001    eq 5, p5
     519 
     520          if (usth_param .eq. "gal") then
     521            Salt_us   =   (log(2.868) - log(1 + SaltMo)) * rCd10n/0.085
     522            Salt_us   = Salt_us * Fac_Mo
     523C +...  Salt_us   :  Extension of  Guyomarc'h & Merindol 1998 with
     524C +...              de Montmollin (1978). Gallee et al. 2001
     525          endif
     526 
     527          if (usth_param .eq. "lis") then !Liston et al. 2007
     528            if(ro__SV(ikl,isn)>300.) then
     529              Salt_us   = 0.005*exp(0.013*ro__SV(ikl,isn))
     530            else
     531              Salt_us   = 0.01*exp(0.003*ro__SV(ikl,isn))
     532            endif
     533          endif
     534 
     535          SnowOK   =  1 -min(1,iabs(isn-isnoSV(ikl))) !Switch new vs old snow
     536 
     537          usthSV(ikl) =     SnowOK *   (Salt_us)
     538     .                + (1.-SnowOK)*    usthSV(ikl)
     539 
     540        END DO
     541
     542
     543 
     544!  Feeback between blowing snow turbulent Scale  u* (commented here
     545!  since ustar is an input variable (not in/out) of inlandsis)
     546!  -----------------------------------------------------------------
     547
     548
     549!           VVa_OK      =  max(0.000001,       VVaSBL(ikl))
     550!           sss__N      =  vonkar      *       VVa_OK
     551!           sss__F      = (sqrCm0(ikl) - psim_z + psim_0)
     552!           usuth0      =  sss__N /sss__F                ! u* if NO Blow. Snow
     553 
     554!           sss__G      =  0.27417     * gravit
     555 
     556! !  ______________               _____
     557! !  Newton-Raphson (! Iteration, BEGIN)
     558! !  ~~~~~~~~~~~~~~               ~~~~~
     559!           DO iit=1,nit
     560!           sss__K      =  gravit      * r_Turb * A_Turb *za__SV(ikl)
     561!      .                                     *rCDmSV(ikl)*rCDmSV(ikl)
     562!      .                           /(1.+0.608*QaT_SV(ikl)-qsnoSV(ikl))
     563!           us_127      =  exp(    SblPom *log(us__SV(ikl)))
     564!           us_227      =  us_127         *    us__SV(ikl)
     565!           us_327      =  us_227         *    us__SV(ikl)
     566!           us_427      =  us_327         *    us__SV(ikl)
     567!           us_527      =  us_427         *    us__SV(ikl)
     568 
     569!           us__SV(ikl) =  us__SV(ikl)
     570!      .    - (  us_527     *sss__F     /sss__N
     571!      .      -  us_427
     572!      .      -  us_227     *qsnoSV(ikl)*sss__K
     573!      .      + (us__SV(ikl)*us__SV(ikl)-usthSV(ikl)*usthSV(ikl))/sss__G)
     574!      .     /(  us_427*5.27*sss__F     /sss__N
     575!      .      -  us_327*4.27
     576!      .      -  us_127*2.27*qsnoSV(ikl)*sss__K
     577!      .      +  us__SV(ikl)*2.0                                 /sss__G)
     578 
     579!           us__SV(ikl)= min(us__SV(ikl),usuth0)
     580!           us__SV(ikl)= max(us__SV(ikl),epsi  )
     581!           rCDmSV(ikl)=     us__SV(ikl)/VVa_OK
     582! ! #AE     sss__F     =     vonkar     /rCDmSV(ikl)
     583!           ENDDO
     584 
     585! !  ______________               ___
     586! !  Newton-Raphson (! Iteration, END  )
     587! !  ~~~~~~~~~~~~~~               ~~~
     588 
     589!           us_127      =  exp(    SblPom *log(us__SV(ikl)))
     590!           us_227      =  us_127         *    us__SV(ikl)
     591 
     592! !  Momentum            Turbulent Scale  u*: 0-Limit in case of no Blow. Snow
     593! !  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     594!           dusuth      =  us__SV(ikl) - usthSV(ikl)       ! u* - uth*
     595!           signus      =  max(sign(unun,dusuth),zero)     ! 1 <=> u* - uth* > 0
     596!           us__SV(ikl) =                                  !
     597!      .                   us__SV(ikl)  *signus  +         ! u* (_BS)
     598!      .                   usuth0                          ! u* (nBS)
     599!      .                            *(1.-signus)           !       
     600
     601
     602
     603
     604!  Blowing Snow        Turbulent Scale ss*
     605!  ---------------------------------------
     606 
     607        hSalSV(ikl) = 8.436e-2  * us__SV(ikl)**SblPom
     608 
     609        if (qsalt_param .eq. "pom") then
     610          qSalSV(ikl) = (us__SV(ikl)**2 - usthSV(ikl)**2) *signus
     611     .               / (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25)
     612        endif
     613 
     614        if (qsalt_param .eq. "bin") then
     615          qSalSV(ikl) = (us__SV(ikl) * us__SV(ikl)
     616     .                -usthSV(ikl) * usthSV(ikl))*signus
     617     .                * 0.535 / (hSalSV(ikl) * gravit)
     618        endif
     619 
     620        qSalSV(ikl) = qSalSV(ikl)/rht_SV(ikl) ! conversion kg/m3 to kg/kg
     621 
     622        ssstar      = rCDmSV(ikl) * (qsnoSV(ikl) - qSalSV(ikl))
     623     .              * r_Turb !Bintanja 2000, BLM
     624!r_Turb compensates for an overestim. of the blown snow part. fall velocity
     625 
     626        uss_SV(ikl) = min(zero    , us__SV(ikl) *ssstar)
     627        uss_SV(ikl) = max(-0.0001 , uss_SV(ikl))   
     628
     629
     630
     631
     632        ENDIF   ! BloMod
     633 
     634C + ------------------------------------------------------
    391635C +--Buffer Layer
    392 C +  ------------
     636C +  -----------------------------------------------------
    393637 
    394638          DO ikl=1,knonv
    … …  
    414658c #NP.         104. *sqrt( max( VV10SV(ikl)-6.0,0.0)))  ! Kotlyakov (1961)
    415659 
    416             density_kotlyakov = .true.
    417 c #AC       density_kotlyakov = .false.  !C.Agosta snow densisty as if BS is on b
     660!          C.Agosta option for snow density, same as for BS i.e.
     661!          is_ok_density_kotlyakov=.false.
    418662c #BS       density_kotlyakov = .false.  !C.Amory BS 2018
    419663C + ...     Fallen Snow Density, Adapted for Antarctica
    420             if (density_kotlyakov) then
     664            if (is_ok_density_kotlyakov) then
    421665                tt_tmp = TaT_SV(ikl)-TfSnow
    422666                !vv_tmp = VV10SV(ikl)
    … …  
    452696!    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    453697 
    454 c #BS       Bros_N      = frsno
    455 c #BS       ro_new      = ro__SV(ikl,max(1,isnoSV(ikl)))
    456 c #BS       ro_new      = max(Bros_N,min(roBdSV,ro_new))
    457 c #BS       Fac         = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
    458 c #BS.                     -roBdSV)/(500.-roBdSV))
    459 c #BS       Fac         = max(0.,min(1.,Fac))
    460 c #BS       dsnbSV(ikl) = Fac*dsnbSV(ikl)
    461 c #BS       Bros_N      = Bros_N     * (1.0-dsnbSV(ikl))
    462 c #BS.                  + ro_new     *      dsnbSV(ikl)
    463 
     698         if (BloMod) then
     699         Bros_N      = frsno
     700         ro_new      = ro__SV(ikl,max(1,isnoSV(ikl)))
     701         ro_new      = max(Bros_N,min(roBdSV,ro_new))
     702         Fac         = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
     703     .               -roBdSV)/(500.-roBdSV))
     704         Fac         = max(0.,min(1.,Fac))
     705         dsnbSV(ikl) = Fac*dsnbSV(ikl)
     706         Bros_N      = Bros_N     * (1.0-dsnbSV(ikl))
     707     .               + ro_new     *      dsnbSV(ikl)
     708         endif
    464709
    465710 
    … …  
    480725     .               max(Spher1*VV__SV(ikl)+Spher2,     !     Sphericity
    481726     .                   Spher3                   ))    !
     727! EV: now control buf_sph_pol and bug_siz_pol in physiq.def
    482728            Buf_G1      = (1. - Polair) *   Buf_G1      ! Temperate Snow
    483      .                        + Polair  *   G1_dSV      ! Polar     Snow
     729     .                        + Polair  *   buf_sph_pol ! Polar Snow
    484730            Buf_G2      = (1. - Polair) *   Buf_G2      ! Temperate Snow
    485      .                        + Polair  *   ADSdSV      ! Polar    Snow
     731     .                        + Polair  *   buf_siz_pol ! Polar Snow
    486732                G1      =                   Buf_G1      ! NO  Blown Snow
    487733                G2      =                   Buf_G2      ! NO  Blown Snow
    488734
    489  
     735
     736
     737            IF (BloMod) THEN
     738
    490739!     S.1. Meme  Type  de Neige  / same Grain Type
    491740!          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    492 c #BS       SameOK  =  max(zero,
    493 c #BS.                     sign(unun,    Buf_G1             *G1_dSV
    494 c #BS.                                 - eps_21                    ))
    495 c #BS       G1same  = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV)
    496 c #BS       G2same  = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV)
     741
     742           SameOK  =  max(zero,
     743     .         sign(unun,    Buf_G1             *G1_dSV
     744     .                            - eps_21                    ))
     745           G1same  = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV)
     746           G2same  = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV)
    497747!           Blowing Snow Properties:                         G1_dSV, ADSdSV
    498748 
    499749!     S.2. Types differents / differents Types
    500750!          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    501 c #BS       typ__1  =  max(zero,sign(unun,epsi-Buf_G1))   ! =1.=> Dendritic
    502 c #BS       zroNEW  =     typ__1  *(1.0-dsnbSV(ikl))      ! fract.Dendr.Lay.
    503 c #BS.              + (1.-typ__1) *     dsnbSV(ikl)       !
    504 c #BS       G1_NEW  =     typ__1  *Buf_G1                 ! G1 of Dendr.Lay.
    505 c #BS.              + (1.-typ__1) *G1_dSV                 !
    506 c #BS       G2_NEW  =     typ__1  *Buf_G2                 ! G2 of Dendr.Lay.
    507 c #BS.              + (1.-typ__1) *ADSdSV                 !
    508 c #BS       zroOLD  = (1.-typ__1) *(1.0-dsnbSV(ikl))      ! fract.Spher.Lay.
    509 c #BS.              +     typ__1  *     dsnbSV(ikl)       !
    510 c #BS       G1_OLD  = (1.-typ__1) *Buf_G1                 ! G1 of Spher.Lay.
    511 c #BS.              +     typ__1  *G1_dSV                 !
    512 c #BS       G2_OLD  = (1.-typ__1) *Buf_G2                 ! G2 of Spher.Lay.
    513 c #BS.              +     typ__1  *ADSdSV                 !
    514 c #BS       SizNEW  =    -G1_NEW  *DDcdSV/G1_dSV          ! Size  Dendr.Lay.
    515 c #BS.               +(1.+G1_NEW         /G1_dSV)         !
    516 c #BS.                  *(G2_NEW  *DScdSV/G1_dSV          !
    517 c #BS.               +(1.-G2_NEW         /G1_dSV)*DFcdSV) !
    518 c #BS       SphNEW  =     G2_NEW         /G1_dSV          ! Spher.Dendr.Lay.
    519 c #BS       SizOLD  =     G2_OLD                          ! Size  Spher.Lay.
    520 c #BS       SphOLD  =     G1_OLD         /G1_dSV          ! Spher.Spher.Lay.
    521 c #BS       Siz_av =     (zroNEW*SizNEW+zroOLD*SizOLD)    ! Averaged Size
    522 c #BS       Sph_av = min( zroNEW*SphNEW+zroOLD*SphOLD     !
    523 c #BS.                   ,  unun)                         ! Averaged Sphericity
    524 c #BS       Den_av = min((Siz_av -(    Sph_av *DScdSV     !
    525 c #BS.                            +(1.-Sph_av)*DFcdSV))   !
    526 c #BS.                 / (DDcdSV -(    Sph_av *DScdSV     !
    527 c #BS.                            +(1.-Sph_av)*DFcdSV))   !
    528 c #BS.                   ,  unun)                         !
    529 c #BS       DendOK  = max(zero,                           !
    530 c #BS.                    sign(unun,     Sph_av *DScdSV   ! Small   Grains
    531 c #BS.                              +(1.-Sph_av)*DFcdSV   ! Faceted Grains
    532 c #BS.                              -    Siz_av        )) !
     751           typ__1  =  max(zero,sign(unun,epsi-Buf_G1))   ! =1.=> Dendritic
     752           zroNEW  =     typ__1  *(1.0-dsnbSV(ikl))      ! fract.Dendr.Lay.
     753     .            + (1.-typ__1) *     dsnbSV(ikl)       !
     754           G1_NEW  =     typ__1  *Buf_G1                 ! G1 of Dendr.Lay.
     755     .            + (1.-typ__1) *G1_dSV                 !
     756           G2_NEW  =     typ__1  *Buf_G2                 ! G2 of Dendr.Lay.
     757     .            + (1.-typ__1) *ADSdSV                 !
     758           zroOLD  = (1.-typ__1) *(1.0-dsnbSV(ikl))      ! fract.Spher.Lay.
     759     .            +     typ__1  *     dsnbSV(ikl)       !
     760           G1_OLD  = (1.-typ__1) *Buf_G1                 ! G1 of Spher.Lay.
     761     .            +     typ__1  *G1_dSV                 !
     762           G2_OLD  = (1.-typ__1) *Buf_G2                 ! G2 of Spher.Lay.
     763     .            +     typ__1  *ADSdSV                 !
     764           SizNEW  =    -G1_NEW  *DDcdSV/G1_dSV          ! Size  Dendr.Lay.
     765     .            +(1.+G1_NEW         /G1_dSV)          !
     766     .                  *(G2_NEW  *DScdSV/G1_dSV        !
     767     .            +(1.-G2_NEW         /G1_dSV)*DFcdSV) !
     768           SphNEW  =     G2_NEW         /G1_dSV          ! Spher.Dendr.Lay.
     769           SizOLD  =     G2_OLD                          ! Size  Spher.Lay.
     770           SphOLD  =     G1_OLD         /G1_dSV          ! Spher.Spher.Lay.
     771           Siz_av  =     (zroNEW*SizNEW+zroOLD*SizOLD)   ! Averaged Size
     772           Sph_av  = min( zroNEW*SphNEW+zroOLD*SphOLD    !
     773     .                 ,  unun)                         ! Averaged Sphericity
     774           Den_av  = min((Siz_av -(    Sph_av *DScdSV    !
     775     .            +(1.-Sph_av)*DFcdSV))                 !
     776     .            / (DDcdSV -(    Sph_av *DScdSV        !
     777     .            +(1.-Sph_av)*DFcdSV))                 !
     778     .                   ,  unun)                       !
     779           DendOK  = max(zero,                           !
     780     .                    sign(unun,     Sph_av *DScdSV   ! Small   Grains
     781     .                              +(1.-Sph_av)*DFcdSV   ! Faceted Grains
     782     .                              -    Siz_av        )) !
    533783C +...      REMARQUE: le  type moyen (dendritique ou non) depend
    534784C +         ^^^^^^^^  de la  comparaison avec le diametre optique
    … …  
    538788C +                   of a recent snow    having zero dendricity
    539789 
    540 c #BS       G1diff  =(   -DendOK *Den_av
    541 c #BS.               +(1.-DendOK)*Sph_av) *G1_dSV
    542 c #BS       G2diff  =     DendOK *Sph_av  *G1_dSV
    543 c #BS.               +(1.-DendOK)*Siz_av
    544 c #BS       G1      =     SameOK *G1same
    545 c #BS.               +(1.-SameOK)*G1diff
    546 c #BS       G2      =     SameOK *G2same
    547 c #BS.               +(1.-SameOK)*G2diff
    548  
     790           G1diff  =(   -DendOK *Den_av
     791     .            +(1.-DendOK)*Sph_av) *G1_dSV
     792           G2diff  =     DendOK *Sph_av  *G1_dSV
     793     .            +(1.-DendOK)*Siz_av
     794           G1      =     SameOK *G1same
     795     .            +(1.-SameOK)*G1diff
     796           G2      =     SameOK *G2same
     797     .            +(1.-SameOK)*G2diff
     798           ENDIF
     799
    549800
    550801 
    … …  
    634885     .                            /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m]
    635886 
    636 
    637887 
    638888          END DO
    … …  
    640890
    641891
    642 ! Snow Pack Discretization
    643 ! ========================
    644 
    645 !            **********
     892! Snow Pack Discretization(option XF in MAR)
     893! ==========================================
     894
     895         
     896      if (discret_xf.AND.klonv.eq.1) then
     897
     898       if(isnoSV(1).ge.1.or.NLaysv(1).ge.1) then
     899C +          **********
     900         call SISVAT_zSn
     901C +          **********
     902       endif
     903      else
     904C +          **********
    646905        call SISVAT_zSn
    647 !            **********
    648 
    649 !            **********
     906C +          **********
     907      endif
     908 
     909C +          **********
    650910! #ve   call SISVAT_wEq('_zSn  ',0)
    651 !            **********
    652 
    653 
     911C +          **********
    654912
    655913! Add a new Snow Layer
    … …  
    664922            TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl))
    665923     .                  + min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl)
    666 
    667924            ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl))
    668925     .                      + Brossv(ikl)     *    NLaysv(ikl)
    … …  
    699956
    700957
    701       END IF
     958      END IF  ! SnoMod
    702959
    703960
    … …  
    740997! =============================
    741998!Etienne: as in inlandis we do not call vgopt, we need to define
    742 !the albedo  alb_SV and to calculate the
     999!the albedo alb_SV and to calculate the
    7431000!absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv
    7441001
    … …  
    8101067
    8111068
    812 ! Aerodynamic Resistance
    813 ! ^^^^^^^^^^^^^^^^^^^^^^
    814 
    815 
    816        DO ikl=1,knonv
    817           ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
    818           rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
    819         END DO
     1069! Aerodynamic Resistance (calculated from drags given by LMDZ)
     1070! Commented because already calculated in surf_inlandsis_mod
     1071! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     1072!       DO ikl=1,knonv
     1073!          ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
     1074!          rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
     1075!        END DO
    8201076
    8211077
    … …  
    8251081
    8261082
    827       if (iflag_tsurf_inlandsis .eq. 0) then
     1083      if (iflag_temp_inlandsis .eq. 0) then
    8281084
    8291085       call SISVAT_TSo
    8301086
    8311087      else
     1088        DO ikl=1,knonv
     1089        Tsf_SV(ikl)=Tsrfsv(ikl)
     1090        END DO
    8321091
    8331092       call SISVAT_TS2
    … …  
    9381197! Surface Temperature
    9391198! ^^^^^^^^^^^^^^^^^^^^
    940 !           Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
    941 
     1199
     1200          IF (iflag_tsurf_inlandsis .EQ. 0) THEN   
     1201
     1202            Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
     1203
     1204          ELSE IF (iflag_tsurf_inlandsis .GT. 0) THEN
    9421205! Etienne: extrapolation from the two uppermost levels:
    9431206
    … …  
    9591222
    9601223
    961         END DO
    962 
     1224         ELSE !(default)
     1225
     1226           Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
     1227
     1228         END IF
     1229
     1230
     1231         END DO
    9631232
    9641233! Snow Pack Properties (sphericity, dendricity, size)
    … …  
    9671236      IF (SnoMod)                                                 THEN
    9681237
    969 !            **********
     1238      if (discret_xf .AND. klonv.eq.1) then
     1239      if(isnoSV(1).ge.1) then
     1240C +          **********
     1241      call SISVAT_GSn
     1242C +          **********
     1243      endif
     1244      else
     1245C +          **********
    9701246        call SISVAT_GSn
    971 !            **********
    972 
    973 !            **********
    974 ! #ve   call SISVAT_wEq('_GSn  ',0)
    975 !            **********
    976 
     1247C +          **********
     1248      endif
    9771249
    9781250
    … …  
    9901262C +--Roughness Length for Momentum
    9911263C +  -----------------------------
     1264
     1265! ETIENNE WARNING: changes have been made wrt original SISVAT
    9921266 
    9931267C +--Land+Sea-Ice / Ice-free Sea Mask
    9941268C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    995         DO ikl=1,klonv
     1269        DO ikl=1,knonv
    9961270          IcIndx(ikl) = 0
    9971271        ENDDO
    9981272        DO isn=1,nsno
    999         DO ikl=1,klonv
     1273        DO ikl=1,knonv
     1274
    10001275          IcIndx(ikl) = max(IcIndx(ikl),
    1001      .                      isn*max(0,
    1002      .                              sign(1,
    1003      .                                   int(ro__SV(ikl,isn)-900.))))
     1276     .                  isn*max(0,
     1277     .                  sign(1,
     1278     .                  int(ro__SV(ikl,isn)-900.))))
    10041279        ENDDO
    10051280        ENDDO
    10061281 
    1007         DO ikl=1,klonv
     1282        DO ikl=1,knonv
    10081283          LISmsk    =     1. ! in inlandsis, land only
    10091284          IceMsk    =     max(0,sign(1   ,IcIndx(ikl)-1)  )
    10101285          SnoMsk    = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
    10111286
    1012  
    1013 
    1014           Z0mLnd      =max( Z0_ICE    ,    5.e-5  )  ! Min set := Z0 on *
    10151287
    10161288C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
    10171289C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
    10181290          Z0m_nu =       5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
    1019  
     1291
    10201292C +--Z0 Saltat.Regime over Snow (Gallee  et al., 2001, BLM 99 (19) p.11)
    10211293C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
     1294
    10221295          u2star =       us__SV(ikl) *us__SV(ikl)
    10231296          Z0mBSn =       u2star      *0.536e-3   -  61.8e-6
    10241297          Z0mBSn =   max(Z0mBS0      ,Z0mBSn)
    1025  
     1298
    10261299C +--Z0 Smooth + Saltat. Regime
    10271300C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
    10281301          Z0enSV(ikl) =  Z0m_nu
    10291302     .                +  Z0mBSn
    1030  
    1031 C +--Rough   Snow Surface Roughness Length (Typical Value)
    1032 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    1033 c #tz     Z0m_Sn =    0.250e-3 ! Andreas 1995, CRREL Report 95-16, fig.1&p.2
    1034                                ! z0r~(10-d)*exp(-vonkar/sqrt(1.5e-03))-5.e-5
    1035           Z0m_Sn =    2.000e-3 ! Calibration    of MAR
    1036 c #TZ     Z0m_Sn =    1.000e-3 ! Exemple Tuning in RACMO
    1037 c #TZ     Z0m_Sn =    0.500e-3 ! Exemple Tuning in MAR
    1038  
     1303
     1304       
     1305! Calculation of snow roughness length
     1306!=====================================
     1307          IF (iflag_z0m_snow .EQ. 0) THEN
     1308
     1309          Z0m_Sn=prescribed_z0m_snow
     1310
     1311          ELSE IF (iflag_z0m_snow .EQ. 1) THEN
     1312
     1313          Z0m_Sn=Z0enSV(ikl)
     1314
     1315          ELSE IF (iflag_z0m_snow .EQ. 2) THEN                             
     1316
    10391317C +--Rough   Snow Surface Roughness Length (Variable Sastrugi Height)
    10401318C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    … …  
    10451323! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    10461324! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013
     1325
     1326         
    10471327          coefa = 0.1658 !0.1862 !Ant
    10481328          coefb = -50.3869 !-55.7718 !Ant
    … …  
    10551335          coefc = log(z03/z02)/(ta3-ta2)
    10561336          coefd = log(z03)-coefc*ta3
     1337
    10571338          if (TaT_SV(ikl) .lt. ta1) then
    10581339            Z0_obs = z01
    … …  
    10661347          endif
    10671348 
    1068 
    1069 ! pour le moment, on choisit une valeur fixe
    1070           Z0_obs = 1.000e-3
    1071  
    1072 cCA       Snow roughness lenght deduced from observations
    1073 cCA       (parametrization if no Blowing Snow)
    1074 cCA       ----------------------------------- C. Agosta 09-2016 -----
    1075 cCA       Substract Z0enSV(ikl) because re-added later in Z0mnSV(ikl)
    1076           Z0m_Sn = Z0_obs - Z0enSV(ikl)
    1077 cCA       -----------------------------------------------------------
    1078  
    1079           param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING
    1080  
     1349          Z0m_Sn=Z0_obs
     1350
     1351
     1352          ELSE
     1353
     1354          Z0m_Sn=0.500e-3  ! default=0.500e-3m (tuning of MAR)
     1355
     1356          ENDIF
     1357 
     1358
     1359
     1360!          param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING
    10811361c #SZ     Z0Sa_N =                   (us__SV(ikl) -0.2)*param   ! 0.0001=TUNING
    10821362c #SZ.           * max(zero,sign(unun,TfSnow-eps9
    … …  
    11091389c #ZN     Z0enSV(ikl) =  max(Z0enSV(ikl), Z0m_nu)
    11101390 
     1391
    11111392C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
    11121393C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^  ams.confex.com/ams/pdfpapers/68601.pdf)
    … …  
    11321413c #ZA     Z0m_Sn =           DDs_SV(ikl)* Z0m_90 / 45.
    11331414c #ZA.         - DDs_SV(ikl)*DDs_SV(ikl)* Z0m_90 /(90.*90.)
    1134  
    1135 C +--Z0  (Erosion)    over Snow (instantaneous or time average)
     1415
     1416
     1417
     1418
     1419C +--Z0  (Erosion)    over Snow (instantaneous)
    11361420C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
    11371421          Z0e_SV(ikl) =  Z0enSV(ikl)
    1138           Z0e_SV(ikl) =  Z0emSV(ikl)
    1139  
    1140 C +--Momentum  Roughness Length
    1141 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^                              ! Contribution of
    1142           Z0mnSV(ikl) =  Z0mLnd                              ! land Form
    1143      .                + (Z0m_Sn                              ! Sastrugi   Form
    1144      .                +  Z0enSV(ikl))   *SnoMsk              ! Snow    Erosion
     1422 
     1423C +--Momentum  Roughness Length (Etienne: changes wrt original SISVAT)
     1424C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                             
     1425          Z0mnSV(ikl) =  Z0m_nu *(1-SnoMsk)                     ! Ice z0
     1426     .                + (Z0m_Sn)*SnoMsk                         ! Snow Sastrugi Form and Snow Erosion
    11451427 
    11461428
    … …  
    11541436c #GL.                     /(920.00                 -600.))) !
    11551437 
    1156 C +--Mom. Roughness Length, Instantaneous OR Box Moving Average in Time
    1157 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     1438C +--Mom. Roughness Length, Instantaneous
     1439C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    11581440          Z0m_SV(ikl) =  Z0mnSV(ikl)                         ! Z0mnSV  instant.
    1159 !          Z0m_SV(ikl) =  Z0mmSV(ikl)                         ! Z0mnSV  Average
    1160  
    1161 C +--Corrected Threshold Friction Velocity before Erosion    ! Marticorena and
    1162 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^    ! Bergametti 1995
    1163 ! not used anymore since Marticorena and Bergametti disabled !CK 18/07/2018
    1164 cc #BS     Z0e_SV(ikl) =   min(Z0m_SV(ikl),Z0e_SV(ikl))       !
    1165 cc #MB     f_eff=    log(0.35*(0.1        /Z0e_SV(ikl))**0.8) ! JGR 100
    1166 cc #MB     f_eff=1.-(log(      Z0m_SV(ikl)/Z0e_SV(ikl)      ))! (20) p. 16420
    1167 cc #MB.            /(max(      f_eff      ,epsi             ))! p.16426 2nd ?
    1168 cc #MB     f_eff=    max(      f_eff      ,epsi              )! CONTROL
    1169 cc #MB     f_eff=1.0   -(1.0 - f_eff)     /5.00               ! TUNING
    1170 cc #MB     f_eff=    min(      f_eff      ,1.00              )!
    1171 cc #MB    usthSV(ikl) =       usthSV(ikl)/f_eff              !
    1172  
    11731441 
    11741442 
    … …  
    11771445 
    11781446          Z0hnSV(ikl) =     Z0mnSV(ikl)/  7.4
    1179 c #SH     Z0hnSV(ikl) =     Z0mnSV(ikl)/100.0
    1180 C +                         Z0h = Z0m  /100.0   over the Sahel
    1181 C +                                            (Taylor & Clark, QJRMS 127,p864)
    1182  
    1183 c #RN     rstar       =     Z0mnSV(ikl) * us__SV(ikl) / akmol
    1184 c #RN     rstar       = max(epsi,min(rstar,thous))
    1185 c #RN     alors       =          log(rstar)
    1186 c #RN     rstar0      = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar))
    1187 c #RN.                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
    1188 c #RN.                *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar))
    1189 c #RN.                + 0.317e0
    1190 c #RN.                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
    1191 c #RN     rstar1      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
    1192 c #RN.                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
    1193 c #RN.                *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar))
    1194 c #RN.                - 0.565
    1195 c #RN.                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
    1196 c #RN     rstar2      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
    1197 c #RN.                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
    1198 c #RN.                *(0.      * max(zero,sign(unun,2.500e0 - rstar))
    1199 c #RN.                - 0.183
    1200 c #RN.                *(unun    - max(zero,sign(unun,2.500e0 - rstar))))
    1201  
    1202 cXF    #RN does not work over bare ice
    1203 cXF    MAR is then too warm and not enough melt
    1204  
    1205 c #RN     if(ro__SV(ikl,isnoSV(ikl))>50
    1206 c #RN.  .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
    1207  
    1208 c #RN     Z0hnSV(ikl) = max(zero
    1209 c #RN.                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))
    1210 c #RN.                * exp(rstar0+rstar1*alors+rstar2*alors*alors)
    1211 c #RN.                * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero
    1212 c #RN.                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
    1213  
    1214 c #RN     endif
     1447 
     1448          IF (is_ok_z0h_rn) THEN
     1449
     1450          rstar       =     Z0mnSV(ikl) * us__SV(ikl) / akmol
     1451          rstar       = max(epsi,min(rstar,R_1000))
     1452          alors       =          log(rstar)
     1453          rstar0      = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar))
     1454     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
     1455     .                *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar))
     1456     .                + 0.317e0
     1457     .                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
     1458          rstar1      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
     1459     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
     1460     .                *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar))
     1461     .                - 0.565
     1462     .                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
     1463          rstar2      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
     1464     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
     1465     .                *(0.      * max(zero,sign(unun,2.500e0 - rstar))
     1466     .                - 0.183
     1467     .                *(unun    - max(zero,sign(unun,2.500e0 - rstar))))
     1468 
     1469         
     1470
     1471!XF    #RN (is_ok_z0h_rn) does not work well over bare ice
     1472!XF    MAR is then too warm and not enough melt
     1473 
     1474         if(ro__SV(ikl,isnoSV(ikl))>50
     1475     .  .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
     1476 
     1477             Z0hnSV(ikl) = max(zero
     1478     .                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))
     1479     .                * exp(rstar0+rstar1*alors+rstar2*alors*alors)
     1480     .                * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero
     1481     .                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
     1482 
     1483          endif
     1484
     1485
     1486          ENDIF
    12151487 
    12161488          Z0h_SV(ikl) =     Z0hnSV(ikl)
    1217 !          Z0h_SV(ikl) =     Z0hmSV(ikl)
    12181489 
    12191490

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/sisvat_bsn.F

    r3792 r4013  
    99C |                                                                        |
    1010C |   SISVAT_bsn computes the snow erosion mass according to both the      |
    11 C |   theoretical maximum erosion amount computed in SISVATesbl and the    |
     11C |   theoretical maximum erosion amount computed in inlandsis and the     |
    1212C |   availability of snow (currently in the uppermost snow layer only)    |
    1313C |                                                                        |
    14 C |   Preprocessing  Option: SISVAT IO (not always a standard preprocess.) |
    15 C |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^                                     |
    16 C |   FILE                 |      CONTENT                                  |
    17 C |   ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
    18 C | # stdout               | #sb: OUTPUT of Snow Erosion                   |
    19 C |                        |      unit  6, SubRoutine  SISVAT_BSn **ONLY** |
    2014C +------------------------------------------------------------------------+
    2115 

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/sisvat_qsn.F

    r3792 r4013  
    6161      use VARxSV
    6262      use VARySV
     63      use surface_data, only: is_ok_slush,opt_runoff_ac
     64
    6365
    6466      IMPLICIT NONE
    … …  
    235237 
    236238      DO ikl=1,knonv
    237        DO isn=min(nsno,isnoSV(ikl)+1),1,-1
     239
     240      DO isn=min(nsno,isnoSV(ikl)+1),1,-1
    238241! EV          DO isn=nsno,1,-1
    239242C +--Energy, store Previous Content
    … …  
    243246     .                + ro__SV(ikl,isn) * Cn_dSV * dTSnow
    244247     .                                           * dzsnSV(ikl,isn)
    245 
    246           Tsave       = TsisSV(ikl,isn)
    247 
    248248          TsisSV(ikl,isn) =                        TfSnow
    249249 
    … …  
    312312          rdzNEW      = WaFrez + rdzsno
    313313          ro__SV(ikl,isn) =      rdzNEW /max(epsi, dzsnSV(ikl,isn))
    314 
    315 ! EV: condition on Enfrez
    316 !          if (EnFrez .eq. 0.) then
    317          
    318           TsisSV(ikl,isn) = Tsave
    319 !          else
    320314          TsisSV(ikl,isn) =      TfSnow
    321315     .                + EnFrez /(Cn_dSV *max(epsi, rdzNEW)        )
    322 !          end if
    323316          EExcsv(ikl) =          EExcsv(ikl)     - EnFrez
    324317          wer_SV(ikl) = WaFrez
    … …  
    499492          rusnew      = rusnSV(ikl) * SWf_SV(ikl)
    500493 
    501           if(isnoSV(ikl)<=1) rusnew = 0.
     494          if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0.
    502495          !if(ivgtSV(ikl)>=1) rusnew = 0.
    503496 
    504497c #EU                        rusnew = 0.
    505 c #AC                        rusnew = 0.
     498c #AC               rusnew = 0.
     499
    506500          RnofSV(ikl) = RnofSV(ikl)
    507501     .                +(rusnSV(ikl) - rusnew     ) / dt__SV
    … …  
    545539        ENDDO
    546540 
    547 C +--Slush Formation (CAUTION: ADD RunOff Possibility before Activation)
     541C +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
    548542C +  ---------------  ^^^^^^^  ^^^
    549543 
    550  
    551 c #SU DO  ikl=1,knonv
    552 c #SU  DO isn=1,isnoSV(ikl)
    553 c #SU     kSlush = min(1,max(0,isn+1-ispiSV(ikl)))        ! Slush Switch
     544      IF (is_ok_slush) THEN
     545
     546      DO  ikl=1,knonv
     547       DO isn=1,isnoSV(ikl)
     548          kSlush = min(1,max(0,isn+1-ispiSV(ikl)))        ! Slush Switch
    554549 
    555550C +--Available Additional Pore   Volume [-]
    556551C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    557 c #SU     PorVol = 1. - ro__SV(ikl,isn)                    ! [--]
    558 c #SU.           *(1. - eta_SV(ikl,isn))/ ro_Ice           !
    559 c #SU.           -      eta_SV(ikl,isn)                    !
    560 c #SU.                 *ro__SV(ikl,isn) / ro_Wat           !
    561 c #SU     PorVol =  max(PorVol          , zero  )          !
    562 c #SU     zWater =      dzsnSV(ikl,isn) * PorVol * 1000.   ! [mm] OR [kg/m2]
    563 c #SU.           * (1. -SWS_SV(ikl)                        ! 0 <=> freezing
    564 c #SU.                *(1 -min(1,iabs(isn-isnoSV(ikl)))))  ! 1 <=> isn=isnoSV
    565 c #SU     zSlush =  min(rusnSV(ikl)     , zWater)          ! [mm] OR [kg/m2]
    566 c #SU     ro_new      =(dzsnSV(ikl,isn) * ro__SV(ikl,isn)  !
    567 c #SU.                 +zSlush                           ) !
    568 c #SU.            / max(dzsnSV(ikl,isn) , epsi           ) !
    569 c #SU     if(ro_new<ro_Ice+20) then ! MAX 940kg/m3         !
    570 c #SU      rusnSV(ikl)  = rusnSV(ikl)          - zSlush    ! [mm] OR [kg/m2]
    571 c #SU      RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)
    572 c #SU      eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn)      !
    573 c #SU.                     *(1.     - eta_SV(ikl,isn)))    !
    574 c #SU.                / max (ro_new , epsi            )    !
    575 c #SU      ro__SV(ikl,isn) =      ro_new                   !
    576 c #SU     endif
    577 c #SU   END DO
    578 c #SU END DO
    579  
     552          PorVol = 1. - ro__SV(ikl,isn)                    ! [--]
     553     .           *(1. - eta_SV(ikl,isn))/ ro_Ice           !
     554     .           -      eta_SV(ikl,isn)                    !
     555     .                 *ro__SV(ikl,isn) / ro_Wat           !
     556          PorVol =  max(PorVol          , zero  )          !
     557          zWater =      dzsnSV(ikl,isn) * PorVol * 1000.   ! [mm] OR [kg/m2]
     558     .           * (1. -SWS_SV(ikl)                        ! 0 <=> freezing
     559     .                *(1 -min(1,iabs(isn-isnoSV(ikl)))))  ! 1 <=> isn=isnoSV
     560          zSlush =  min(rusnSV(ikl)     , zWater)          ! [mm] OR [kg/m2]
     561          ro_new      =(dzsnSV(ikl,isn) * ro__SV(ikl,isn)  !
     562     .                 +zSlush                           ) !
     563     .            / max(dzsnSV(ikl,isn) , epsi           ) !
     564          if(ro_new<ro_Ice+20) then ! MAX 940kg/m3         !
     565           rusnSV(ikl)  = rusnSV(ikl)          - zSlush    ! [mm] OR [kg/m2]
     566           RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)
     567           eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn)      !
     568     .                     *(1.     - eta_SV(ikl,isn)))    !
     569     .                / max (ro_new , epsi            )    !
     570           ro__SV(ikl,isn) =      ro_new                   !
     571          endif
     572        END DO
     573      END DO
     574      END IF
    580575 
    581576C +--Impact of the Sublimation/Deposition on the Surface Mass Balance

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/sisvat_tso.F

    r3792 r4013  
    132132
    133133      integer   nt_srf,it_srf,itEuBk          ! HL: Surface Scheme
    134       parameter(nt_srf=10)                    !
     134      parameter(nt_srf=10)                     ! 10 before
    135135      real      agpsrf,xgpsrf,dt_srf,dt_ver   !
    136136      real      etaBAK(knonv)                 !
    … …  
    153153C +                                           ! including   Snow Melt  Energy
    154154 
    155  
     155C +-- Initilialisation of local arrays
     156C +   ================================
     157        DO ikl=1,knonv
     158
     159          mu_sno(ikl)=0.
     160          mu__dz(ikl,:)=0.   
     161          dtC_sv(ikl,:)=0.     
     162          IRs__D(ikl)=0.                 
     163          dIRsdT(ikl)=0.                 
     164          f_HSHL(ikl)=0.                 
     165          dRidTs(ikl)=0.                 
     166          HS___D(ikl)=0.               
     167          f___HL(ikl)=0.             
     168          HL___D(ikl)=0.           
     169          TSurf0(ikl)=0.     
     170          qsatsg(ikl)=0.
     171          dqs_dT(ikl)=0.                 
     172          Psi(ikl)=0.               
     173          RHuSol(ikl)=0.                 
     174          Diag_A(ikl,:)=0.     
     175          Diag_B(ikl,:)=0.     
     176          Diag_C(ikl,:)=0.     
     177          Term_D(ikl,:)=0.     
     178          Aux__P(ikl,:)=0.     
     179          Aux__Q(ikl,:)=0.     
     180          etaBAK(ikl)=0.               
     181          etaNEW(ikl)=0.               
     182          etEuBk(ikl)=0.                 
     183          fac_dt(ikl)=0.
     184          faceta(ikl)=0. 
     185          PsiArg(ikl)=0.
     186          SHuSol(ikl)=0. 
     187
     188        END DO
     189
     190       
    156191
    157192C +--Heat Conduction Coefficient (zero in the Layers over the highest one)
    … …  
    336371C +--Snow highest Layer (dummy!)
    337372C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
    338         isl=  min(isnoSV(1)+1,nsno)
    339         DO ikl=1,knonv
     373
     374        !EV!isl=  min(isnoSV(1)+1,nsno)
     375
     376        DO ikl=1,knonv
     377! EV try to calculate isl at the ikl grid point
     378          isl=  min(isnoSV(ikl)+1,nsno)
     379
    340380          Elem_A          =  dtC_sv(ikl,isl)  *mu__dz(ikl,isl)
    341381          Elem_C          =  0.
    … …  
    384424c    .                                   / den_qs              !
    385425c         qsatsg(ikl) = .0038 *        exp(arg_qs)             !
    386 
    387426!          sp = (pst_SV(ikl) + ptopSV) * 10.
    388427
    389           sp=ps__SV(ikl)
     428          !sp=ps__SV(ikl)
     429          ! Etienne: in the formula herebelow sp should be in hPa, not
     430          ! in Pa so I divide by 100.
     431          sp=ps__SV(ikl)/100.
    390432          psat_ice = 6.1070 * exp(6150. *(1./273.16 -
    391433     .                                              1./TsisSV(ikl,isl)))
    … …  
    399441            qsatsg(ikl) = 0.622 * psat_wat / (sp - 0.378 * psat_wat)
    400442          endif
     443          QsT_SV(ikl)=qsatsg(ikl)
    401444
    402445c         dqs_dT(ikl) = qsatsg(ikl)* 4099.2   /(den_qs *den_qs)!
    403446          fac_dt(ikl) = f_HSHL(ikl)/(ro_Wat   * dz_dSV(0))     !
    404447        END DO
     448
     449
    405450 
    406451C +--Surface: Latent    Heat Flux: Surface    Relative Humidity
    … …  
    410455     .                             /(   1.0-xgpsrf**nt_srf)    !
    411456              dt_srf       = agpsrf                            !
    412               dt_ver       = 0.                                !
     457              dt_ver       = 0.               
     458
    413459            DO ikl=1,knonv
    414               isl          =          isnoSV(ikl)              !
     460              isl          =          isnoSV(ikl)             
     461              ist          = max(0,isotSV(ikl)-100*isnoSV(ikl))! 0 if    H2O                         
     462              ist__s       = min(1,ist)                                                                             
    415463              etaBAK(ikl)  = max(epsi,eta_SV(ikl ,isl))        !
    416464              etaNEW(ikl)  =          etaBAK(ikl)              !
    417465              etEuBk(ikl)  =          etaNEW(ikl)              !
    418             END DO                                             !
     466            END DO     
     467
     468        if(ist__s==1) then ! to reduce computer time                                                 
     469                                          !
    419470        DO it_srf=1,nt_srf                                     !
    420471              dt_ver       = dt_ver     +dt_srf                !
    … …  
    458509            END DO                                             !
    459510              dt_srf      =      dt_srf         * xgpsrf       !
    460         END DO                                                 !
     511        END DO         
     512
     513
     514        endif                                       !
    461515 
    462516C +--Surface: Latent    Heat Flux: Soil/Water Surface Contributions
    … …  
    579633       
    580634      END DO
     635
     636
    581637 
    582638C +--Temperature Limits (avoids problems in case of no Snow Layers)
    … …  
    584640        DO ikl=     1,knonv
    585641           isl              = isnoSV(ikl)
    586           dTSurf            = TsisSV(ikl,isl) -     TSurf0(ikl)
     642
     643           dTSurf            = TsisSV(ikl,isl) -     TSurf0(ikl)
    587644          TsisSV(ikl,isl)   = TSurf0(ikl) + sign(1.,dTSurf) ! 180.0 dgC/hr
    588645     .              * min(abs(dTSurf),5.e-2*dt__SV)         ! =0.05 dgC/s
    … …  
    602659C +--Update Surface    Fluxes
    603660C +  ========================
    604  
     661       
     662
     663
    605664        DO ikl=      1,knonv
    606665          isl         = isnoSV(ikl)
    … …  
    613672        END DO
    614673
    615  
    616  
    617674      return
    618675      end

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/sisvat_zsn.F

    r3792 r4013  
    5252      use VARxSV
    5353      use VARySV
     54      use surface_data, only: ok_zsn_ii
    5455
    5556      IMPLICIT NONE
    … …  
    716717        END DO
    717718
     719
     720C +--Search new Ice/Snow Interface (option II in MAR)
     721C +  ===============================================
     722
     723        IF (ok_zsn_ii) THEN
     724       
     725        DO ikl=1,knonv
     726          iiceSV(ikl) =  0
     727        END DO
     728 
     729        DO ikl=1,knonv
     730        DO   isn=1,isnoSV(ikl)
     731          OK_ICE      = max(zero,sign(unun,ro__SV(ikl,isn)-ro_ice+20.))
     732     .                * max(zero,sign(unun,dzsnSV(ikl,isn)-epsi))
     733          iiceSV(ikl) = (1.-OK_ICE)       *iiceSV(ikl)
     734     .                +     OK_ICE        *isn
     735        END DO
     736        END DO
     737
     738        END IF
    718739 
    719740      return

  • LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/surf_inlandsis_mod.F90

    r3792 r4013  
    11MODULE surf_inlandsis_mod
    22
    3   IMPLICIT NONE
    4 
     3    IMPLICIT NONE
     4
     5CONTAINS
     6
     7
     8SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, &
     9            rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
     10            precip_rain, precip_snow, &
     11            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
     12            rugos, snow_cont_air, alb_soil, alt, slope, cloudf, &
     13            radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
     14            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
     15            runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, &
     16            tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf)
     17
     18        ! |                                                                        |
     19        ! |   SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow|
     20        ! |                              (INLANDSIS)                               |
     21        ! |   SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)         |
     22        ! |   surface scheme of the Modele Atmospherique Regional (MAR)            |
     23        ! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
     24        ! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
     25        ! |   Updated by Etienne Vignon, Cecile Agosta                             |
     26        ! |                                                                        |
     27        ! +------------------------------------------------------------------------+
     28        ! |
     29        ! |   In the current setup, SISVAT is used only to model the land ice      |
     30        ! |   part of the surface; hence it is called with the compressed variables|
     31        ! |   from pbl_surface, and only by the surf_landice routine.              |
     32        ! |                                                                        |
     33        ! |   In this interface it is assumed that the partitioning of the soil,   |
     34        ! |   and hence the number of grid points is constant during a simulation, |
     35        ! |   hence eg. snow properties remain stored in the global SISVAT         |
     36        ! |   variables between the calls and don't need to be handed over as      |
     37        ! |   arguments. When the partitioning is supposed to change, make sure to |
     38        ! |   update the variables.                                                |
     39        ! |                                                                        |
     40        ! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV ...)                    |
     41        ! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |
     42        ! |                                                                        |
     43        ! +------------------------------------------------------------------------+
     44
     45        USE dimphy
     46        USE VAR_SV
     47        USE VARdSV
     48        USE VARxSV
     49        USE VARySV
     50        USE VARtSV
     51        USE VARphy
     52        USE surface_data, only : iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor
     53
     54        IMPLICIT NONE
     55
     56        ! +--INTERFACE Variables
     57        ! +  ===================
     58        !    include  "dimsoil.h"
     59
     60        ! +--Global Variables
     61        ! +  ================
     62        ! Input Variables for SISVAT
     63        INTEGER, INTENT(IN) :: knon
     64        INTEGER, INTENT(IN) :: itime
     65        REAL, INTENT(IN) :: dtime
     66        LOGICAL, INTENT(IN) :: debut     ! true if first step
     67        LOGICAL, INTENT(IN) :: lafin     ! true if last step
     68
     69        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i     ! Index Decompression
     70        REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
     71        REAL, DIMENSION(klon), INTENT(IN) :: rmu0      ! cos sol. zenith angle
     72        REAL, DIMENSION(klon), INTENT(IN) :: swdown    !
     73        REAL, DIMENSION(klon), INTENT(IN) :: lwdown    !
     74        REAL, DIMENSION(klon), INTENT(IN) :: albedo_old
     75        REAL, DIMENSION(klon), INTENT(IN) :: pexner    ! Exner potential
     76        REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
     77        REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo
     78        REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay
     79        REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf
     80        REAL, DIMENSION(klon), INTENT(IN) :: rugos
     81        REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air
     82        REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope
     83        REAL, DIMENSION(klon), INTENT(IN) :: alt       ! surface elevation
     84        REAL, DIMENSION(klon), INTENT(IN) :: cloudf
     85        REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
     86        REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
     87        REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh
     88        REAL, DIMENSION(klon), INTENT(IN) :: ustar   ! friction velocity
     89
     90        ! Variables exchanged between LMDZ and SISVAT
     91        REAL, DIMENSION(klon), INTENT(IN) :: radsol    ! Surface absorbed rad.
     92        REAL, DIMENSION(klon), INTENT(INOUT) :: snow      ! Tot snow mass [kg/m2]
     93        REAL, DIMENSION(klon), INTENT(INOUT) :: zfra      ! snwo surface fraction [0-1]
     94        REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
     95        REAL, DIMENSION(klon), INTENT(OUT) :: qsol      ! Soil Water Content
     96        REAL, DIMENSION(klon), INTENT(INOUT) :: z0m    ! Momentum Roughn Lgt
     97        REAL, DIMENSION(klon), INTENT(INOUT) :: z0h    ! Momentum Roughn Lgt
     98
     99        ! Output Variables for LMDZ
     100        REAL, DIMENSION(klon), INTENT(OUT) :: alb1      ! Albedo SW
     101        REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW
     102        REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo
     103        REAL, DIMENSION(klon), INTENT(OUT) :: emis_new  ! Surface Emissivity
     104        REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff
     105        REAL, DIMENSION(klon), INTENT(OUT) :: ffonte    ! enthalpy flux due to surface melting
     106        REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte   ! water flux due to surface melting
     107        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s   ! d/dT sens. ht flux
     108        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l   ! d/dT latent ht flux
     109        REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens  ! Sensible ht flux
     110        REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat   ! Latent heat flux
     111        REAL, DIMENSION(klon), INTENT(OUT) :: evap      ! Evaporation
     112        REAL, DIMENSION(klon), INTENT(OUT) :: erod      ! Erosion of surface snow (flux)
     113        REAL, DIMENSION(klon), INTENT(OUT) :: agesno    ! Snow age (top layer)
     114        REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature
     115        REAL, DIMENSION(klon), INTENT(OUT) :: qsurf     ! Surface Humidity
     116
     117        ! Specific INLANDIS outputs
     118        REAL, DIMENSION(klon), INTENT(OUT) :: qsnow     ! Total H2O snow[kg/m2]
     119        REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt   ! Snow height (m)
     120        REAL, DIMENSION(klon), INTENT(OUT) :: to_ice    ! Snow passed to ice
     121        REAL, DIMENSION(klon), INTENT(OUT) :: sissnow   ! Snow in model (kg/m2)
     122
     123        ! +--Internal  Variables
     124        ! +  ===================
     125
     126        CHARACTER(len = 20) :: fn_outfor ! Name for output file
     127        CHARACTER (len = 80)              :: abort_message
     128        CHARACTER (len = 20)              :: modname = 'surf_inlandsis_mod'
     129
     130        INTEGER :: i, ig, ikl, isl, isn, nt
     131        INTEGER :: gp_outfor, un_outfor
     132        REAL, PARAMETER :: f1 = 0.5
     133        REAL, PARAMETER :: sn_upp = 10000., sn_low = 500.
     134        REAL, PARAMETER :: sn_add = 400., sn_div = 2.
     135        ! snow mass upper,lower limit,
     136        ! added mass/division lowest layer
     137        REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6
     138        REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3
     139        ! Parameters for drainage
     140        ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
     141        ! +...        Run Off Parameters
     142        ! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)
     143        ! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)
     144
     145        REAL, DIMENSION(klon) :: eps0SL          ! surface Emissivity
     146        REAL :: zsigma, Ua_min, Us_min, lati
     147        REAL, PARAMETER :: cdmax=0.05
     148        REAL :: lambda          ! Par. soil discret.
     149        REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2         ! Soil layer thicknesses
     150        !$OMP THREADPRIVATE(dz1,dz2)
     151        LOGICAL, SAVE :: firstcall
     152        !$OMP THREADPRIVATE(firstcall)
     153
     154        INTEGER :: iso
     155        LOGICAL :: file_exists
     156        CHARACTER(len = 20) :: fichnom
     157        LOGICAL :: is_init_domec
     158        ! CA initialization
     159        ! dz_profil_15 : 1 m in 15 layers [m]
     160        real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, &
     161                                                0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/)
     162        ! mean_temp : mean annual surface temperature [K]
     163        real, dimension(klon) :: mean_temp
     164        ! mean_dens : mean surface density [kg/m3]
     165        real, dimension(klon) :: mean_dens
     166        ! lat_scale : temperature lapse rate against latitude [K degree-1]
     167        real :: lat_scale
     168        ! sh_scale : temperature lapse rate against altitude [K km-1]
     169        real :: sh_scale
     170        ! variables for density profile
     171        ! E0, E1 : exponent
     172        real :: E0, E1
     173        ! depth at which 550 kg m-3 is reached [m]
     174        real :: z550
     175        ! depths of snow layers
     176        real :: depth, snow_depth, distup
     177        ! number of initial snow layers
     178        integer :: nb_snow_layer
     179        ! For density calc.
     180        real :: alpha0, alpha1, ln_smb
     181        ! theoritical densities [kg m-3]
     182        real :: rho0, rho1, rho1_550
     183        ! constants for density profile
     184        ! C0, C1 : constant, 0.07 for z <= 550 kg m-3
     185        real, parameter :: C0 = 0.07
     186        real, parameter :: C1 = 0.03
     187        ! rho_i : ice density [kg m-3]
     188        real, parameter :: rho_ice = 917.
     189        ! E_c : activation energy [J mol-1]
     190        real, parameter :: E_c = 60000.
     191        ! E_g : activation energy [J mol-1]
     192        real, parameter :: E_g = 42400.
     193        ! R : gas constant [J mol-1 K-1]
     194        real, parameter :: R = 8.3144621
     195
     196     
     197     
     198
     199
     200        ! + PROGRAM START
     201        ! + -----------------------------------------
     202
     203        zsigma = 1000.
     204        dt__SV = dtime
     205
     206        IF (debut) THEN
     207            firstcall = .TRUE.
     208            INI_SV = .false.
     209        ELSE
     210            firstcall = .false.
     211            INI_SV = .true.
     212        END IF
     213
     214        IF (ok_outfor) THEN
     215            un_outfor = 51    ! unit number for point output file
     216            gp_outfor = 1    ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC
     217            fn_outfor = 'outfor_SV.dat'
     218        END IF
     219
     220        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     221        ! + INITIALISATION: BEGIN +++
     222        ! + -----------------------------------------
     223        IF (firstcall) THEN
     224
     225            ! +--Array size
     226            ! +  -----------------------
     227
     228            klonv = klon
     229            knonv = knon
     230                write(*, *) 'ikl, lon and lat in INLANDSIS'
     231
     232            DO ikl = 1, knon
     233                i=ikl2i(ikl)
     234                write(*, *) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i)
     235            END DO
     236
     237            ! +--Variables initizialisation
     238            ! +  ---------------------------
     239
     240            CALL INIT_VARtSV
     241            CALL INIT_VARxSV
     242            CALL INIT_VARySV
     243
     244
     245
     246            ! +--Surface Fall Line Slope
     247            ! +  -----------------------
     248            IF (SnoMod)  THEN
     249                DO ikl = 1, knon
     250                    slopSV(ikl) = slope(ikl)
     251                    SWf_SV(ikl) = &   ! Normalized Decay of the
     252                            exp(-dt__SV             &   ! Surficial Water Content
     253                                    / (c1_zuo                &   !(Zuo and Oerlemans 1996,
     254                                            + c2_zuo * exp(-c3_zuo * abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
     255                END DO
     256            END IF
     257
     258
     259
     260            ! +--Soil layer thickness . Compute soil discretization (as for LMDZ)
     261            ! +  ----------------------------------------------------------------
     262            !        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
     263            CALL get_soil_levels(dz1, dz2, lambda)
     264
     265            lambSV = lambda
     266            dz1_SV(1:knon, 1:) = 0.
     267            dz2_SV(1:knon, 1:) = 0.
     268
     269            DO isl = -nsol, 0
     270                dz_dSV(isl) = 0.5e-3 * dz2(1 - isl)           ! Soil layer thickness
     271                DO ikl = 1, knon
     272                    dz1_SV(ikl, isl) = dz1(1 - isl)    !1.e-3*
     273                    dz2_SV(ikl, isl) = dz2(1 - isl)    !1.e-3*
     274                END DO
     275            END DO
     276
     277
     278            ! Set variables
     279            ! =============
     280            DO ikl = 1, knon
     281                ! LSmask : Land/Sea Mask
     282                LSmask(ikl) = 1
     283                ! isotSV : Soil Type -> 12 = ice
     284                isotSV(ikl) = 12
     285                ! iWaFSV : Soil Drainage (1,0)=(y,n)
     286                iWaFSV(ikl) = 1
     287                ! eps0SL : Surface Emissivity
     288                eps0SL(ikl) = 1.
     289                ! alb0SV : Soil Albedo
     290                alb0SV(ikl) = alb_soil(ikl)
     291                ! Tsf_SV : Surface Temperature, must be bellow freezing
     292                Tsf_SV(ikl) = min(temp_air(ikl), TfSnow)
     293            END DO
     294
     295            ! +--Initialization of soil and snow variables in case startsis is not read
     296            ! +  ----------------------------------------------------------------------
     297
     298            is_init_domec=.FALSE.
     299
     300
     301            IF (is_init_domec) THEN
     302            ! Coarse initilization inspired from vertcical profiles at Dome C,
     303            ! Antarctic Plateaui (10m of snow, 19 levels)
     304
     305                 DO ikl = 1,knon
     306! + Soil
     307                 DO isl =   -nsol,0   
     308                   TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2)   !temp_air(ikl)
     309                   !tsoil(ikl,1-isl)   Soil Temperature
     310                   !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
     311                   eta_SV(ikl,isl) = epsi           !etasoil(ikl,1-isl) Soil Water[m3/m3]
     312                   ro__SV(ikl,isl) = rhoIce         !rosoil(ikl,1-isl) volumic mass
     313                 END DO   
     314
     315
     316           ! Snow
     317                 isnoSV(ikl) = 19
     318                 istoSV(ikl, 1:isnoSV(ikl)) = 100
     319                 ro__SV(ikl, 1:isnoSV(ikl)) = 350.
     320                 eta_SV(ikl, 1:isnoSV(ikl)) = epsi
     321                 TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2)
     322                 G1snSV(ikl, 1:isnoSV(ikl)) = 99.
     323                 G2snSV(ikl, 1:isnoSV(ikl)) = 2.
     324                 agsnSV(ikl, 1:isnoSV(ikl)) = 50.
     325                 dzsnSV(ikl, 19) = 0.015
     326                 dzsnSV(ikl, 18) = 0.015
     327                 dzsnSV(ikl, 17) = 0.020
     328                 dzsnSV(ikl, 16) = 0.030
     329                 dzsnSV(ikl, 15) = 0.040
     330                 dzsnSV(ikl, 14) = 0.060
     331                 dzsnSV(ikl, 13) = 0.080
     332                 dzsnSV(ikl, 12) = 0.110
     333                 dzsnSV(ikl, 11) = 0.150
     334                 dzsnSV(ikl, 10) = 0.200
     335                 dzsnSV(ikl, 9) = 0.300
     336                 dzsnSV(ikl, 8) = 0.420
     337                 dzsnSV(ikl, 7) = 0.780
     338                 dzsnSV(ikl, 6) = 1.020
     339                 dzsnSV(ikl, 5) = 0.980
     340                 dzsnSV(ikl, 4) = 1.020
     341                 dzsnSV(ikl, 3) = 3.970
     342                 dzsnSV(ikl, 2) = 1.020
     343                 dzsnSV(ikl, 1) = 1.020
     344
     345                 END DO
     346            ELSE
     347
     348            ! Initilialisation with climatological temperature and density
     349            ! profiles as in MAR. Methodology developed by Cecile Agosta
     350 
     351            ! initialize with 0., for unused snow layers
     352            dzsnSV = 0.
     353            G1snSV = 0.
     354            G2snSV = 0.
     355            istoSV = 0
     356            TsisSV = 0.
     357
     358
     359            ! initialize mean variables (unrealistic)
     360            mean_temp = TfSnow
     361            mean_dens = 300.
     362            ! loop on grid cells
     363            DO ikl = 1, knon
     364                lati=rlat(ikl2i(ikl))
     365                ! approximations for mean_temp and mean_dens
     366                ! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1)
     367                ! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1
     368                ! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica,
     369                ! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed.
     370                if (lati > 60.) then
     371                    ! CA todo : add longitude bounds
     372                    ! Greenland mean temperature : function of altitude and latitude
     373                    ! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1
     374                    lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9
     375                    lat_scale = max(min(lat_scale, 0.9), 0.75)
     376                    ! sh_scale equals the environmental lapse rate : 6.5 °C km-1
     377                    sh_scale = 6.5
     378                    mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.)
     379                    ! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051
     380                    mean_dens(ikl) = 315.
     381                else if (lati < -60.) then
     382                    ! Antarctica mean temperature : function of altitude and latitude
     383                    ! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1
     384                    lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3
     385                    lat_scale = max(min(lat_scale, 1.3), 0.6)
     386                    ! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1
     387                    sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5
     388                    sh_scale = max(min(sh_scale, 9.8), 6.5)
     389                    mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.)
     390                    ! Antarctica surface density : function of mean annual temperature
     391                    ! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013)
     392                    ! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m.
     393                    ! Weinhart et al 2020  https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201
     394                    ! 320 kg m-3 is reached at Dome A, 4100 m a.s.l.
     395                    ! Dome C : st_ant_param(3233, -75.1) = -47.7
     396                    ! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7
     397                    mean_dens(ikl) =  (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450.
     398                    mean_dens(ikl) = min(450., max(320., mean_dens(ikl)))
     399                else
     400
     401                !    write(*, *) 'Attention: temperature initialization is only defined for Greenland and Antarctica'
     402
     403                     mean_dens(ikl) =350.
     404                     mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2)
     405
     406                !abort_message='temperature initialization is only defined for Greenland and Antarctica'
     407                !CALL abort_physic(modname,abort_message,1)
     408
     409                end if
    5410 
    6   CONTAINS
    7 
    8                                        
    9 
    10   SUBROUTINE surf_inlandsis(knon,rlon,rlat, ikl2i, itime, dtime, debut, lafin, &
    11              rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
    12              precip_rain, precip_snow, precip_snow_adv, snow_adv, &
    13              zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
    14              rugos, snow_cont_air, alb_soil, slope, cloudf, &
    15              radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
    16              AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
    17              runoff_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &     
    18              tsurf_new, alb1, alb2, alb3, &
    19              emis_new, z0m, z0h, qsurf)       
    20                                                                              
    21 ! +------------------------------------------------------------------------+   
    22 ! |                                                                        |   
    23 ! |   SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow|
    24 ! |                              (INLANDSIS)                               |
    25 ! |   SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)         |   
    26 ! |   surface scheme of the Modele Atmospherique Regional (MAR)            |
    27 ! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
    28 ! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
    29 ! |           Update Etienne Vignon, LMD, Novembre 2020                    |
    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 ! |             SnoMod: Snow Pack is set up when .T.                       | 
    46 ! |             reaLBC: Update Bound.Condit.when .T.                       |   
    47 ! |                                                                        | 
    48 ! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV)                        | 
    49 ! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |   
    50 ! |                                                                        | 
    51 ! |   Preprocessing  Option: SISVAT PHYSICS                                |   
    52 ! |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^^^^^^                                |   
    53 ! | #                       #HY                                            |   
    54 ! | #                       #SN: Snow         Model                        |   
    55 ! | #                       #BS: Blowing Snow Parameterization             |   
    56 ! +------------------------------------------------------------------------+
    57        
    58     USE dimphy                                           
    59     USE VAR_SV
    60     USE VARdSV
    61     USE VARxSV
    62     USE VARySV
    63     USE VARtSV
    64     USE VARphy
    65     USE surface_data, only: iflag_tsurf_inlandsis,SnoMod,BloMod,ok_outfor                                       
    66 
    67     IMPLICIT NONE   
    68                                                                
    69 ! +--INTERFACE Variables                                                     
    70 ! +  ===================
    71                                            
    72 !    include  "dimsoil.h"                                     
    73                                                                    
    74 
    75 ! +--Global Variables                                                         
    76 ! +  ================ 
    77 ! Input Variables for SISVAT
    78     INTEGER,               INTENT(IN)      :: knon
    79     INTEGER,               INTENT(IN)      :: itime   
    80     REAL,                  INTENT(IN)      :: dtime
    81     LOGICAL,               INTENT(IN)      :: debut     ! true if first step
    82     LOGICAL,               INTENT(IN)      :: lafin     ! true if last step
    83 
    84     INTEGER, DIMENSION(klon), INTENT(IN)   :: ikl2i     ! Index Decompression
    85     REAL, DIMENSION(klon), INTENT(IN)      :: rlon, rlat
    86     REAL, DIMENSION(klon), INTENT(IN)      :: rmu0      ! cos sol. zenith angle
    87     REAL, DIMENSION(klon), INTENT(IN)      :: swdown    !
    88     REAL, DIMENSION(klon), INTENT(IN)      :: lwdown    !
    89     REAL, DIMENSION(klon), INTENT(IN)      :: albedo_old   
    90     REAL, DIMENSION(klon), INTENT(IN)      :: pexner    ! Exner potential
    91     REAL, DIMENSION(klon), INTENT(IN)      :: precip_rain, precip_snow
    92     REAL, DIMENSION(klon), INTENT(IN)      :: precip_snow_adv, snow_adv
    93                                                         !Snow Drift
    94     REAL, DIMENSION(klon), INTENT(IN)      :: zsl_height, wind_velo
    95     REAL, DIMENSION(klon), INTENT(IN)      :: temp_air, spechum, ps,p1lay
    96     REAL, DIMENSION(klon), INTENT(IN)      :: dens_air, tsurf           
    97     REAL, DIMENSION(klon), INTENT(IN)      :: rugos,snow_cont_air
    98     REAL, DIMENSION(klon), INTENT(IN)      :: alb_soil, slope
    99     REAL, DIMENSION(klon), INTENT(IN)      :: cloudf   
    100     REAL, DIMENSION(klon), INTENT(IN)      :: AcoefH, AcoefQ
    101     REAL, DIMENSION(klon), INTENT(IN)      :: BcoefH, BcoefQ
    102     REAL, DIMENSION(klon), INTENT(IN)      :: cdragm, cdragh
    103     REAL, DIMENSION(klon), INTENT(IN)      :: ustar   ! friction velocity
    104 
    105 ! Variables exchanged between LMDZ and SISVAT
    106     REAL, DIMENSION(klon), INTENT(IN)      :: radsol    ! Surface absorbed rad.
    107     REAL, DIMENSION(klon), INTENT(INOUT)   :: snow      ! Tot snow mass [kg/m2]
    108     REAL, DIMENSION(klon), INTENT(INOUT)   :: zfra      ! snwo surface fraction [0-1]
    109     REAL, DIMENSION(klon,nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
    110     REAL, DIMENSION(klon), INTENT(OUT)       :: qsol      ! Soil Water Content 
    111     REAL, DIMENSION(klon), INTENT(INOUT)     :: z0m    ! Momentum Roughn Lgt
    112     REAL, DIMENSION(klon), INTENT(INOUT)     :: z0h    ! Momentum Roughn Lgt
    113 
    114 
    115 ! Output Variables for LMDZ
    116     REAL, DIMENSION(klon), INTENT(OUT)     :: alb1      ! Albedo SW
    117     REAL, DIMENSION(klon), INTENT(OUT)     :: alb2,alb3 ! Albedo NIR and LW
    118     REAL, DIMENSION(klon), INTENT(OUT)     :: emis_new  ! Surface Emissivity
    119     REAL, DIMENSION(klon), INTENT(OUT)     :: runoff_lic ! Runoff           
    120     REAL, DIMENSION(klon), INTENT(OUT)     :: dflux_s   ! d/dT sens. ht flux 
    121     REAL, DIMENSION(klon), INTENT(OUT)     :: dflux_l   ! d/dT latent ht flux
    122     REAL, DIMENSION(klon), INTENT(OUT)     :: fluxsens  ! Sensible ht flux   
    123     REAL, DIMENSION(klon), INTENT(OUT)     :: fluxlat   ! Latent heat flux
    124     REAL, DIMENSION(klon), INTENT(OUT)     :: evap      ! Evaporation
    125     REAL, DIMENSION(klon), INTENT(OUT)     :: agesno    ! Snow age (top layer)
    126     REAL, DIMENSION(klon), INTENT(OUT)     :: tsurf_new ! Surface Temperature
    127     REAL, DIMENSION(klon), INTENT(OUT)     :: qsurf     ! Surface Humidity
    128 
    129 ! Specific INLANDIS outputs
    130 
    131     REAL, DIMENSION(klon), INTENT(OUT)     :: qsnow     ! Total H2O snow[kg/m2]
    132     REAL, DIMENSION(klon), INTENT(OUT)     :: snowhgt   ! Snow height (m)
    133     REAL, DIMENSION(klon), INTENT(OUT)     :: to_ice    ! Snow passed to ice
    134     REAL, DIMENSION(klon), INTENT(OUT)     :: sissnow   ! Snow in model (kg/m2)
    135 
    136                                                                          
    137 
    138                                                                    
    139 ! +--Internal  Variables                                                       
    140 ! +  ===================                                         
    141 
    142     CHARACTER(len=20)               :: fn_outfor ! Name for output file
    143     INTEGER                         :: i, ig, ikl, isl, isn, nt
    144     INTEGER                         :: gp_outfor, un_outfor
    145     REAL, PARAMETER                 :: f1=0.5
    146     REAL, PARAMETER                 :: sn_upp=5000.,sn_low=500.
    147     REAL, PARAMETER                 :: sn_add=400.,sn_div=2.
    148                                              ! snow mass upper,lower limit,
    149                                              ! added mass/division lowest layer
    150     REAL, PARAMETER                 :: c1_zuo=12.960e+4, c2_zuo=2.160e+6
    151     REAL, PARAMETER                 :: c3_zuo=1.400e+2,  czemin=1.e-3 
    152                                              ! Parameters for drainage
    153 ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
    154 ! +...        Run Off Parameters                                             
    155 ! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)   
    156 ! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)             
    157 
    158     REAL, DIMENSION(klon)           :: eps0SL          ! surface Emissivity
    159     REAL                            :: zsigma, Ua_min, Us_min
    160     REAL                            :: lambda          ! Par. soil discret.
    161     REAL, DIMENSION(nsoilmx), SAVE  :: dz1,dz2         ! Soil layer thicknesses
    162 !$OMP THREADPRIVATE(dz1,dz2)
    163     LOGICAL, SAVE                   :: firstcall
    164 !$OMP THREADPRIVATE(firstcall)               
    165 
    166                
    167                                        
    168 ! +--Internal Variables
    169 ! +  ==================
    170 
    171     INTEGER                         ::  iso
    172     LOGICAL                         ::  file_exists
    173     CHARACTER(len=20)               ::  fichnom
    174 !========================================================================
    175 
    176       PRINT*, 'je rentre dans inlandsis'
    177 
    178       zsigma=1000.
    179       dt__SV=dtime
    180      
    181    
    182 
    183 !     write(*,*)'Start of simulation? ',debut        !hj
    184 
    185       IF (debut) THEN
    186         firstcall=.TRUE.
    187         INI_SV=.false.
    188 
    189       ELSE
    190         firstcall=.false.
    191         INI_SV=.true.
    192       END IF
    193 
    194 
    195 
    196 
    197        IF (ok_outfor) THEN
    198         un_outfor=51                 ! unit number for point output file
    199         gp_outfor= 1        ! grid point number for point output
    200         fn_outfor='outfor_SV.dat'
    201        END IF
    202 
    203 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    204 
    205 
    206 
    207 
    208 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    209 ! + INITIALISATION: BEGIN +++
    210 ! + -------------------------
    211 ! +
    212 ! + Compute soil discretization (as for LMDZ)
    213 ! + ----------------------------------------- 
    214       IF (firstcall) THEN
    215 
    216 ! +--Array size
    217      klonv=klon
    218      knonv=knon
    219 
    220 
    221         write(*,*)'klon',klon,'klonv',klonv,'knon',knon,'nsol',nsol,'nsno',nsno
    222 
    223 
    224         CALL INIT_VARtSV
    225         CALL INIT_VARxSV
    226         CALL INIT_VARySV
    227 
    228         eps0SL(:)=0.
    229        
    230        
    231 ! +--Soil layer thickness                                                   
    232 ! +  ----------------------- 
    233 !        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
    234         CALL get_soil_levels(dz1,dz2,lambda)
    235 
    236        
    237         lambSV=lambda
    238         dz1_SV(1:knon,1:) = 0.     
    239         dz2_SV(1:knon,1:) = 0.
    240                      
    241         DO isl =   -nsol,0   
    242           dz_dSV(isl) = 0.5e-3*dz2(1-isl)           ! Soil layer thickness
    243           DO ikl=1,knon
    244             dz1_SV(ikl,isl) = dz1(1-isl)    !1.e-3*
    245             dz2_SV(ikl,isl) = dz2(1-isl)    !1.e-3*
    246           END DO
     411                ! mean_temp is defined for ice ground only
     412                mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2)
     413
     414                ! Soil layers
     415                ! ===========
     416                DO isl = -nsol, 0
     417                    ! TsisSV : Temperature [K]
     418                    TsisSV(ikl, isl) = mean_temp(ikl)
     419                    ! eta_SV : Soil Water [m3/m3]
     420                    eta_SV(ikl, isl) = epsi
     421                    ! ro__SV : Volumic Mass [kg/m3]
     422                    ro__SV(ikl, isl) = rhoIce
     423                END DO
     424
     425                ! Snow layers
     426                ! ===========
     427                ! snow_depth : initial snow depth
     428                snow_depth = 20.
     429                ! nb_snow_layer : initial nb of snow layers
     430                nb_snow_layer = 15
     431                ! isnoSV : total nb of snow layers
     432                isnoSV(ikl) = nb_snow_layer
     433                ! depth : depth of each layer
     434                depth = snow_depth
     435                do isl = 1, nb_snow_layer
     436                    ! dzsnSV : snow layer thickness
     437                    dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1))
     438                    ! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical
     439                    G1snSV(ikl, isl) = 99.
     440                    ! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size
     441                    G2snSV(ikl, isl) = 3.
     442                    agsnSV(ikl, isl) = 0.
     443                    istoSV(ikl, isl) = 0
     444                    ! eta_SV : Liquid Water Content [m3/m3]
     445                    eta_SV(ikl, isl) = 0.
     446                    ! distance to surface
     447                    depth = depth - dzsnSV(ikl,isl) / 2.
     448                    distup = min(1., max(0., depth / snow_depth))
     449                    ! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom)
     450                    TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup**2) + mean_temp(ikl) * distup**2
     451                    ! firn density : densification formulas from :
     452                    ! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/)
     453                    ! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306)
     454                    ! Integration of the steady state equation
     455                    ! ln_smb approximated as a function of temperature
     456                    ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.)
     457                    ! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1
     458                    alpha0 = max(1.435 - 0.151 * ln_smb, 0.25)
     459                    alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25)
     460                    E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0
     461                    E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1
     462                    z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0
     463                    rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
     464                    rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice
     465                    if (depth <= z550) then
     466                        ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
     467                    else
     468                        ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice
     469                    end if
     470                    depth = depth - dzsnSV(ikl,isl) / 2.
     471                   
     472                end do
     473
     474            END DO
     475
     476            END IF
     477
     478
     479            ! + Numerics paramaters, SISVAT_ini
     480            ! +  ----------------------
     481            CALL SISVAT_ini(knon)
     482
     483
     484            ! +--Read restart file
     485            ! +  =================================================
     486
     487            INQUIRE(FILE = "startsis.nc", EXIST = file_exists)
     488            IF (file_exists) THEN
     489                CALL sisvatetat0("startsis.nc", ikl2i)
     490            END IF
     491
     492
     493
     494            ! +--Output ascii file
     495            ! +  =================================================
     496
     497            ! open output file
     498            IF (ok_outfor) THEN
     499                open(unit = un_outfor, status = 'replace', file = fn_outfor)
     500                ikl = gp_outfor     ! index sur la grille land ice
     501                write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl))
     502                write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] &
     503                        & G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]'
     504            END IF
     505
     506        END IF  ! firstcall
     507        ! +
     508        ! +  +++  INITIALISATION:  END  +++
     509        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     510
     511
     512
     513        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     514        ! + READ FORCINGS
     515        ! + ------------------------
     516
     517        ! + Update Forcings for SISVAT given by the LMDZ model.
     518        ! +
     519        DO ikl = 1, knon
     520
     521            ! +--Atmospheric Forcing                                    (INPUT)
     522            ! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^
     523            za__SV(ikl) = zsl_height(ikl)               ! surface layer height (fisr model level) [m]
     524            Ua_min = 0.2 * sqrt(za__SV(ikl))            !
     525            VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
     526            TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
     527            ExnrSV(ikl) = pexner(ikl)                   ! Exner potential
     528            rhT_SV(ikl) = dens_air(ikl)                 ! Air density
     529            QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity
     530            ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
     531            p1l_SV(ikl) = p1lay(ikl)                    ! lowest atm. layer press[Pa]
     532
     533            ! +--Surface properties
     534            ! +  ^^^^^^^^^^^^^^^^^^
     535
     536            Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
     537            Z0h_SV(ikl) = z0h(ikl)                      ! Moment.Roughn.L.
     538
     539            ! +--Energy Fluxes                                          (INPUT)
     540            ! +  ^^^^^^^^^^^^^                                           ^^^^^
     541            coszSV(ikl) = max(czemin, rmu0(ikl))         ! cos(zenith.Dist.)
     542            sol_SV(ikl) = swdown(ikl)                   ! downward Solar
     543            IRd_SV(ikl) = lwdown(ikl)                   ! downward IR
     544            rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.
     545
     546            ! +--Water  Fluxes                                          (INPUT)
     547            ! +  ^^^^^^^^^^^^^                                           ^^^^^
     548            drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
     549            dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]
     550
     551            ! #BS    dbs_SV(ikl) = blowSN(i,j,n)
     552            ! dbs_SV = Maximum potential erosion amount [kg/m2]
     553            ! => Upper bound for eroded snow mass
     554            !        uss_SV(ikl) = SLussl(i,j,n) ! u*qs* (only for Tv in sisvatesbl.f)
     555            ! #BS  if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then
     556            ! #BS    dsnbSV(ikl) =1.0-min(qsHY(i,j,kB)     !BS neglib. at kb ~100 magl)
     557            ! #BS.                        /max(qshy(i,j,mz),eps9),unun)
     558            ! #BS    dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl))
     559            ! #BS    dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl)))
     560            ! #BS  else
     561            ! #BS    dsnbSV(ikl) = 0.
     562            ! #BS  endif
     563            !      dsnbSV is the drift fraction of deposited snow updated in sisvat.f
     564            !      will be used for characterizing the Buffer Layer
     565            !      (see update of  Bros_N, G1same, G2same, zroOLD, zroNEW)
     566            ! #BS  if(n==1) qbs_HY(i,j) = dsnbSV(ikl)
     567            qsnoSV(ikl) = snow_cont_air(ikl)
     568
     569
     570
     571            ! +--Soil/BL                                      (INPUT)
     572            ! +  ^^^^^^^                                       ^^^^^
     573            alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo
     574            AcoHSV(ikl) = AcoefH(ikl)
     575            BcoHSV(ikl) = BcoefH(ikl)
     576            AcoQSV(ikl) = AcoefQ(ikl)
     577            BcoQSV(ikl) = BcoefQ(ikl)
     578            cdH_SV(ikl) = min(cdragh(ikl),cdmax)
     579            cdM_SV(ikl) = min(cdragm(ikl),cdmax)
     580            rcdmSV(ikl) = sqrt(cdM_SV(ikl))
     581            Us_min = 0.01
     582            us__SV(ikl) = max(Us_min, ustar(ikl))
     583            ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6))
     584            rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6))
     585
     586            ! +--Energy Fluxes                                          (INPUT/OUTPUT)
     587            ! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^
     588            !IF (.not.firstcall) THEN
     589            Tsrfsv(ikl)  = tsurf(ikl)                     !hj 12 03 2010
     590            cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness
     591            !END IF
     592
     593         END DO
     594
     595        !
     596        ! +  +++  READ FORCINGS:  END  +++
     597        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     598 
     599
     600        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     601        ! +--SISVAT EXECUTION
     602        ! +  ----------------
     603
     604        call  INLANDSIS(SnoMod, BloMod, 1)
     605
     606
     607       
     608        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     609        ! + RETURN RESULTS
     610        ! + --------------
     611        ! + Return (compressed) SISVAT variables to LMDZ
     612        ! +
     613        DO  ikl = 1, knon                  ! use only 1:knon (actual ice sheet..)
     614            dflux_s(ikl) = dSdTSV(ikl)         ! Sens.H.Flux T-Der.
     615            dflux_l(ikl) = dLdTSV(ikl)         ! Latn.H.Flux T-Der.
     616            fluxsens(ikl) = HSs_sv(ikl)         ! HS
     617            fluxlat(ikl) = HLs_sv(ikl)         ! HL
     618            evap(ikl) = -1*HLs_sv(ikl) / LHvH2O  ! Evaporation
     619            erod(ikl) = 0.
     620
     621            IF (BloMod) THEN
     622                ! + Blowing snow
     623
     624                !       SLussl(i,j,n)= 0.
     625                ! #BS   SLussl(i,j,n)=                     !Effective erosion
     626                ! #BS. (- dbs_ER(ikl))/(dt*rhT_SV(ikl))    !~u*qs* from previous time step
     627                ! #BS   blowSN(i,j,n)=  dt*uss_SV(ikl)     !New max. pot. Erosion [kg/m2]
     628                ! #BS.                    *rhT_SV(ikl)     !(further bounded in sisvat_bsn.f)
     629                ! #BS  erprev(i,j,n) =     dbs_Er(ikl)/dt__SV
     630                erod(ikl) = dbs_Er(ikl) / dt__SV
     631            ENDIF
     632
     633            ! +   Check snow thickness,  substract if too thick, add if too thin
     634
     635            sissnow(ikl) = 0.  !()
     636            DO  isn = 1, isnoSV(ikl)
     637                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
     638            END DO
     639
     640            IF (sissnow(ikl) .LE. sn_low) THEN  !add snow
     641                IF (isnoSV(ikl).GE.1) THEN
     642                    dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi)
     643                    toicSV(ikl) = toicSV(ikl) - sn_add
     644                ELSE
     645                    write(*, *) 'Attention, bare ice... point ', ikl
     646                    isnoSV(ikl) = 1
     647                    istoSV(ikl, 1) = 0
     648                    ro__SV(ikl, 1) = 350.
     649                    dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi)  ! 1.
     650                    eta_SV(ikl, 1) = epsi
     651                    TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2)
     652                    G1snSV(ikl, 1) = 0.
     653                    G2snSV(ikl, 1) = 0.3
     654                    agsnSV(ikl, 1) = 10.
     655                    toicSV(ikl) = toicSV(ikl) - sn_add
     656                END IF
     657            END IF
     658
     659            IF (sissnow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
     660                dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div
     661                toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div
     662            END IF
     663
     664            sissnow(ikl) = 0.
     665            qsnow(ikl) = 0.
     666            snow(ikl) = 0.
     667            snowhgt(ikl) = 0.
     668
     669            DO  isn = 1, isnoSV(ikl)
     670                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
     671                snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn)
     672                ! Etienne: check calc qsnow
     673                qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn)
     674            END DO
     675
     676            zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0)
     677            ! Etienne: comment following line
     678            ! snow(ikl)    = sissnow(ikl)+toicSV(ikl)
     679            snow(ikl) = sissnow(ikl)
     680
     681            to_ice(ikl) = toicSV(ikl)
     682            runoff_lic(ikl) = RnofSV(ikl)    ! RunOFF: intensity (flux due to melting + liquid precip)
     683            fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing
     684            ffonte(ikl)=fqfonte(ikl)*Lf_H2O
     685
     686            qsol(ikl) = 0.
     687            DO  isl = -nsol, 0
     688                tsoil(ikl, 1 - isl) = TsisSV(ikl, isl)       ! Soil Temperature
     689                ! Etienne: check calc qsol
     690                qsol(ikl) = qsol(ikl)                      &
     691                        + eta_SV(ikl, isl) * dz_dSV(isl)
     692            END DO
     693            agesno(ikl) = agsnSV(ikl, isnoSV(ikl))        !          [day]
     694
     695            alb1(ikl) = alb1sv(ikl)             ! Albedo VIS
     696!            alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl)                   &
     697!                    & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1
     698            alb2(ikl)=alb2sv(ikl)
     699            ! Albedo NIR
     700            alb3(ikl) = alb3sv(ikl)             ! Albedo FIR
     701            ! 6 band Albedo
     702            alb6(ikl,:)=alb6sv(ikl,:)
     703
     704            tsurf_new(ikl) = Tsrfsv(ikl)
     705
     706            qsurf(ikl) = QsT_SV(ikl)
     707            emis_new(ikl) = eps0SL(ikl)
     708            z0m(ikl) = Z0m_SV(ikl)
     709            z0h(ikl) = Z0h_SV(ikl)
     710
    247711
    248712        END DO
    249713
    250 
    251 
    252 
    253         DO ikl=1,knon     
    254 
    255 
    256 ! Initialise variables
    257          
    258           ispiSV(ikl)             = 0
    259           iiceSV(ikl)             = 0 
    260           rusnSV(ikl)             = 0.   
    261           toicSV(ikl)             = 0.     
    262           isnoSV(ikl)             = 0.       ! # snow layers                           
    263           istoSV(ikl,:)           = 0.
    264           eta_SV(ikl,:)           = 0.     
    265           TsisSV(ikl,:)           = 0.
    266           ro__SV(ikl,:)           = 0.       
    267           G1snSV(ikl,:)           = 0. 
    268           G2snSV(ikl,:)           = 0.
    269           agsnSV(ikl,:)           = 0.
    270           dzsnSV(ikl,:)           = 0.
    271           zzsnsv(ikl,:)           = 0.                                           
    272           BufsSV(ikl)             = 0.   
    273           qsnoSV(ikl)             = 0.     ! BL snow content 
    274           zWEcSV(ikl)             = 0.
    275           dbs_SV(ikl)             = 0.
    276           dsnbSV(ikl)             = 0.
    277           esnbSV(ikl)             = 0.
    278           BrosSV(ikl)             = 0.
    279           BG1sSV(ikl)             = 0.         
    280           BG2sSV(ikl)             = 0.
    281           SWS_SV(ikl)             = 0.
    282           RnofSV(ikl)             = 0.     ! RunOFF Intensity
    283           RRs_SV(ikl)             = 0.
    284           DDs_SV(ikl)             = 0.
    285           VVs_SV(ikl)             = 0.
    286           cld_SV(ikl)             = 0.     
    287           uts_SV(ikl)             = 0.     ! u*T*  arbitrary 
    288           uqs_SV(ikl)             = 0.     ! u*q*    "
    289           uss_SV(ikl)             = 0.     ! u*s*    "
    290           LMO_SV(ikl)             = 0.
    291 
    292 
    293 ! Set variables
    294 
    295           LSmask(ikl) = 1                          ! Land/Sea   Mask   
    296           isotSV(ikl) = 12                         ! Soil       Type  -> 12= ice 
    297           iWaFSV(ikl) = 1                          ! Soil Drainage                                     
    298           eps0SL(ikl )= 1.
    299           alb0SV(ikl) = alb_soil(ikl)                 ! Soil Albedo       
    300           Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
    301           Z0h_SV(ikl) = z0h(ikl)                      ! heat Roughn.L.
    302 
    303 ! + Soil Upward IR Flux, Water Fluxes, roughness length 
    304           IRs_SV(ikl) =                               &
    305               -eps0SL(ikl)* StefBo*(temp_air(ikl)**4)   ! Upward IR Flux   
    306           Tsf_SV(ikl) = min(temp_air(ikl),TfSnow)       
    307          
    308 ! + Soil
    309         DO isl =   -nsol,0   
    310           TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2)   !temp_air(ikl)  !tsoil(ikl,1-isl)   Soil Temperature
    311           !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
    312           eta_SV(ikl,isl) = epsi                        !etasoil(ikl,1-isl) Soil Water[m3/m3]
    313           ro__SV(ikl,isl) = rhoIce                         !rosoil(ikl,1-isl)  volumic mass
    314         END DO     
    315 
    316 
    317 
    318 !! Initialise with snow
    319           !  G1snSV(ikl,0)          = 0.                   !      [-]     
    320           !  G2snSV(ikl,0)          = 1.6                  ! [-] [0.0001 m]
    321           !  dzsnSV(ikl,0)          = dz_dSV(0)            !           [m]
    322 
    323 
    324           ! if (snow(ikl) .GT. 0.) then
    325           !   isnoSV(ikl)             = 1       ! snow layers                           
    326           !   istoSV(ikl,1:nsno)      = 0     ! 0,...,5 :   Snow     History (see istdSV data)
    327           !   eta_SV(ikl,1:nsno)      = epsi     
    328           !   TsisSV(ikl,1:nsno)      = tsoil(ikl,1)         
    329           !   ro__SV(ikl,1:nsno)      = 350.0       
    330           !   G1snSV(ikl,1:nsno)      = 0.   !      [-] 
    331           !   G2snSV(ikl,1:nsno)      = 1.6   !     [-] [0.0001 m]
    332           !   agsnSV(ikl,1:nsno)      = 50.   !          [day]
    333           !   dzsnSV(ikl,1)           = snow(ikl)/max(ro__SV(ikl,1),epsi) ![m]
    334 ! ! ecrete si trop de neige:
    335 !              IF (snow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
    336 !               dzsnSV(ikl,1)      = dzsnSV(ikl,1)/sn_div
    337 !               toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div
    338 !              END IF
    339 !            zzsnsv(ikl,1)      =  dzsnSV(ikl,1)                    ! Total snow pack thickness
    340 !          endif
    341 
    342 
    343 ! Initialise la neige avec un profil de densité prochde des conditions de Dôme C (~10m de neige avec 19 niveaux) (Etienne):
    344           isnoSV(ikl)                    = 19
    345           istoSV(ikl,1:isnoSV(ikl))      = 100
    346           ro__SV(ikl,1:isnoSV(ikl))      = 350.     
    347           eta_SV(ikl,1:isnoSV(ikl))      = epsi
    348           TsisSV(ikl,1:isnoSV(ikl))      = min(tsoil(ikl,1),TfSnow-0.2)
    349           G1snSV(ikl,1:isnoSV(ikl))      = 0
    350           G2snSV(ikl,1:isnoSV(ikl))      = 1.6
    351           agsnSV(ikl,1:isnoSV(ikl))      = 50.   
    352           dzsnSV(ikl,19)                  = 0.015
    353           dzsnSV(ikl,18)                  =0.015
    354           dzsnSV(ikl,17)                  =0.020
    355           dzsnSV(ikl,16)                  =0.030
    356           dzsnSV(ikl,15)                  =0.040
    357           dzsnSV(ikl,14)                  =0.060
    358           dzsnSV(ikl,13)                  =0.080
    359           dzsnSV(ikl,12)                  =0.110
    360           dzsnSV(ikl,11)                  =0.150
    361           dzsnSV(ikl,10)                  =0.200
    362           dzsnSV(ikl,9)                   =0.300
    363           dzsnSV(ikl,8)                   =0.420
    364           dzsnSV(ikl,7)                   =0.780
    365           dzsnSV(ikl,6)                   =1.020
    366           dzsnSV(ikl,5)                   =0.980
    367           dzsnSV(ikl,4)                   =1.020
    368           dzsnSV(ikl,3)                   =3.970
    369           dzsnSV(ikl,2)                   =1.020
    370           dzsnSV(ikl,1)                   =0.100
    371 
    372 
    373         END DO                                           
    374 
    375 ! +--Surface Fall Line Slope                                                   
    376 ! +  ----------------------- 
    377         IF (SnoMod)  THEN               
    378           DO ikl=1,knon 
    379             slopSV(ikl) = slope(ikl)
    380             SWf_SV(ikl) =             &   ! Normalized Decay of the 
    381               exp(-dt__SV             &   ! Surficial Water Content 
    382               /(c1_zuo                &   !(Zuo and Oerlemans 1996, 
    383             +c2_zuo*exp(-c3_zuo*abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
    384           END DO                                     
    385         END IF                           
    386 
    387 ! + SISVAT_ini (as for use with MAR, but not computing soil layers)
    388 ! + -------------------------------------------------------------
    389 !        write(*,'(/a)') 'Start SISVAT initialization: SISVAT_ini'
    390         CALL SISVAT_ini(knon)
    391 
    392 
    393 ! +--Read restart file
    394 ! +  =================================================   
    395        
    396         INQUIRE(FILE="startsis.nc", EXIST=file_exists)
    397         IF (file_exists) THEN
    398         CALL sisvatetat0("startsis.nc",ikl2i)
     714            IF (ok_outfor) THEN
     715             ikl= gp_outfor
     716            write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++'
     717            write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl)
     718            write(un_outfor, *) dzsnSV(ikl, :)
     719            write(un_outfor, *) TsisSV(ikl, :)
     720            write(un_outfor, *) ro__SV(ikl, :)
     721            write(un_outfor, *) eta_SV(ikl, :)
     722            write(un_outfor, *) G1snSV(ikl, :)
     723            write(un_outfor, *) G2snSV(ikl, :)
     724            write(un_outfor, *) agsnSV(ikl, :)
     725            write(un_outfor, *) istoSV(ikl, :)
     726            write(un_outfor, *) DOPsnSV(ikl, :)
     727        ENDIF
     728
     729
     730
     731        ! +  -----------------------------
     732        ! +  END --- RETURN RESULTS
     733        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     734        IF (lafin) THEN
     735            fichnom = "restartsis.nc"
     736            CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat)
     737
     738            IF (ok_outfor) THEN
     739                close(unit = un_outfor)
     740            END IF
    399741        END IF
    400        
    401        
    402        
    403 ! +--Output ascii file
    404 ! +  =================================================   
    405                
    406        
    407        
    408         ! open output file
    409         IF (ok_outfor) THEN
    410           open(unit=un_outfor,status='replace',file=fn_outfor)         
    411           ikl=gp_outfor     ! index sur la grille land ice
    412           write(un_outfor,*) fn_outfor, ikl, dt__SV   
    413           write(un_outfor,*) 'nsnow - albedo - z0m - z0h , dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46] &
    414            & G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]'
    415 
    416         END IF
    417  
    418       END IF  ! firstcall                       
    419 ! +                               
    420 ! +  +++  INITIALISATION:  END  +++                               
    421 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    422 
    423 
    424 
    425 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    426 ! + READ FORCINGS 
    427 ! + ------------------------
    428 
    429 ! + Update Forcings for SISVAT given by the LMDZ model.
    430 ! +
    431       DO ikl=1,knon
    432 
    433 ! +--Atmospheric Forcing                                    (INPUT)           
    434 ! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^
    435         zSBLSV      = 1000.                         ! [m]               
    436         za__SV(ikl) = zsl_height(ikl)               ! surface layer height (fisr model level) [m]
    437         Ua_min      = epsi                          !
    438         Ua_min      = 0.2 * sqrt(za__SV(ikl)   )    !                   
    439         VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
    440         TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
    441         ExnrSV(ikl) = pexner(ikl)                   ! Exner potential         
    442         rhT_SV(ikl) = dens_air(ikl)                 ! Air density           
    443         QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity
    444         ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
    445         p1l_SV(ikl) = p1lay(ikl)                    ! lowest atm. layer press[Pa]
    446 
    447 ! +--Surface properties
    448 ! +  ^^^^^^^^^^^^^^^^^^
    449 
    450         Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
    451         Z0h_SV(ikl) = z0h(ikl)                      ! Moment.Roughn.L.
    452 
    453 ! +--Energy Fluxes                                          (INPUT)           
    454 ! +  ^^^^^^^^^^^^^                                           ^^^^^             
    455         coszSV(ikl) = max(czemin,rmu0(ikl))         ! cos(zenith.Dist.) 
    456         sol_SV(ikl) = swdown(ikl)                   ! downward Solar 
    457         IRd_SV(ikl) = lwdown(ikl)                   ! downward IR   
    458         rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.   
    459 
    460 ! +--Water  Fluxes                                          (INPUT)           
    461 ! +  ^^^^^^^^^^^^^                                           ^^^^^             
    462         drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
    463         dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]
    464 !c #BS  dbsnow      = -SLussl(i,j,n)                ! Erosion   
    465 !c #BS.               *dtPhys     *rhT_SV(ikl) /ro_Wat                   
    466 !c #BS  dsnbSV(ikl) = snow_adv(ikl)  ! min(max(zero,dbsnow)             
    467 !c #BS.                    /    max(epsi,d_snow),unun)                   
    468 !c #BS  dbs_SV(ikl) = snow_cont_air(ikl)
    469 !c #BS                  blowSN(i,j,n)               !          [kg/m2] 
    470                                                                              
    471 ! +--Soil/BL                                      (INPUT)           
    472 ! +  ^^^^^^^                                       ^^^^^           
    473         alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo
    474         AcoHSV(ikl) = AcoefH(ikl) 
    475         BcoHSV(ikl) = BcoefH(ikl)                     
    476         AcoQSV(ikl) = AcoefQ(ikl) 
    477         BcoQSV(ikl) = BcoefQ(ikl)             
    478         cdH_SV(ikl) = cdragh(ikl)     
    479         cdM_SV(ikl) = cdragm(ikl)     
    480         Us_min      = 0.01
    481         us__SV(ikl) = max(Us_min, ustar(ikl))
    482         ram_sv(ikl) = 1./(cdragm(ikl)*max(VV__SV(ikl),eps6))
    483         rah_sv(ikl) = 1./(cdragh(ikl)*max(VV__SV(ikl),eps6))   
    484 
    485 ! +--Energy Fluxes                                          (INPUT/OUTPUT)   
    486 ! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^   
    487         IF (.not.firstcall) THEN 
    488         Tsf_SV(ikl) = tsurf(ikl)                     !hj 12 03 2010
    489         cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness         
    490         END IF
    491  
    492 
    493       END DO
    494 
    495 !                           
    496 ! +  +++  READ FORCINGS:  END  +++   
    497 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    498 
    499 
    500 
    501 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    502 ! +--SISVAT EXECUTION                                                         
    503 ! +  ----------------                                                         
    504 
    505       call  INLANDSIS(SnoMod,BloMod,1)
    506 
    507 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    508 ! + RETURN RESULTS 
    509 ! + --------------
    510 ! + Return (compressed) SISVAT variables to LMDZ             
    511 ! +
    512       DO  ikl=1,knon                  ! use only 1:knon (actual ice sheet..)
    513         runoff_lic(ikl)    = RnofSV(ikl)*dtime   ! RunOFF: intensity* time step
    514         dflux_s(ikl)       = dSdTSV(ikl)         ! Sens.H.Flux T-Der.
    515         dflux_l(ikl)       = dLdTSV(ikl)         ! Latn.H.Flux T-Der.
    516         fluxsens(ikl)      = HSs_sv(ikl)         ! HS                 
    517         fluxlat(ikl)       = HLs_sv(ikl)         ! HL                 
    518         evap(ikl)          = HLs_sv(ikl)/LHvH2O  ! Evaporation 
    519         snow(ikl)          = 0.
    520         snowhgt(ikl)       = 0.
    521         qsnow(ikl)         = 0.
    522         qsol(ikl)          = 0.
    523 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    524 ! +
    525 ! +   Check snow thickness,  substract if too thick   (commended by etienne: add if too thin)
    526 
    527         sissnow(ikl)       = 0.  !()
    528       DO  isn = 1,isnoSV(ikl)                                               
    529         sissnow(ikl)       = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)   
    530       END DO
    531 
    532        IF (sissnow(ikl) .LE. sn_low) THEN  !add snow
    533        IF (isnoSV(ikl).GE.1) THEN
    534          dzsnSV(ikl,1)      = dzsnSV(ikl,1) + sn_add/max(ro__SV(ikl,1),epsi) 
    535          toicSV(ikl)        = toicSV(ikl)   - sn_add
    536 !       ELSE
    537 !         write(*,*) 'Attention, bare ice... point ',ikl
    538 !         isnoSV(ikl)        = 1
    539 !         istoSV(ikl,1)      = 0
    540 !         ro__SV(ikl,1)      = 350.     
    541 !         dzsnSV(ikl,1)      = sn_add/max(ro__SV(ikl,1),epsi)  ! 1.
    542 !         eta_SV(ikl,1)      = epsi
    543 !         TsisSV(ikl,1)      = min(TsisSV(ikl,0),TfSnow-0.2)   
    544 !         G1snSV(ikl,1)      = 0.
    545 !         G2snSV(ikl,1)      = 0.3
    546 !         agsnSV(ikl,1)      = 10.   
    547 !         toicSV(ikl)        = toicSV(ikl)   - sn_add
    548        END IF
    549        END IF
    550 
    551       IF (sissnow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
    552         dzsnSV(ikl,1)      = dzsnSV(ikl,1)/sn_div
    553         toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div
    554       END IF
    555 
    556         sissnow(ikl)       = 0.  !()
    557 
    558         DO  isn = 1,isnoSV(ikl)                                               
    559         sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)           
    560         snowhgt(ikl) = snowhgt(ikl)+dzsnSV(ikl,isn)             
    561         qsnow(ikl)   = qsnow(ikl)+1e03*eta_SV(ikl,isn)*dzsnSV(ikl,isn)   
    562         END DO
    563 
    564         ! Etienne: pourquoi ajouter toicSV ici? Pour bilan d'eau?
    565         snow(ikl)    = sissnow(ikl)+toicSV(ikl)
    566         to_ice(ikl)  = toicSV(ikl)
    567 
    568 
    569       DO  isl =   -nsol,0   
    570         tsoil(ikl,1-isl)   = TsisSV(ikl,isl)       ! Soil Temperature 
    571         qsol(ikl)          = qsol(ikl)                      &   
    572                             +eta_SV(ikl,isl) * dz_dSV(isl) 
    573       END DO                                               
    574         agesno(ikl)        = agsnSV(ikl,isnoSV(ikl))        !          [day]
    575 
    576         alb1(ikl)          = alb1sv(ikl)             ! Albedo VIS 
    577         alb2(ikl)          = ((So1dSV-f1)*alb1sv(ikl)                   &
    578      &                       +So2dSV*alb2sv(ikl)+So3dSV*alb3sv(ikl))/f1   
    579                                                      ! Albedo NIR
    580         alb3(ikl)          = alb3sv(ikl)             ! Albedo FIR 
    581 
    582         tsurf_new(ikl)     =Tsrfsv(ikl)
    583 
    584         zfra(ikl)          = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
    585         qsurf(ikl)         = QaT_SV(ikl)
    586         emis_new(ikl)      = eps0SL(ikl) 
    587         z0m(ikl)           = Z0m_SV(ikl)
    588         z0h(ikl)           = Z0h_SV(ikl)
    589 
    590       END DO  ! ikl
    591      
    592 
    593 
    594 
    595 
    596 
    597 ! write variables in output file
    598 
    599       IF (ok_outfor) THEN
    600         ikl=gp_outfor
    601 
    602 !        write(un_outfor,*) 'nsnow [-,1], dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46]'
    603 !        write(un_outfor,*) 'G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]'
    604         write(un_outfor,*) '+++++++++++++++++++++++++++++++++++++++++++++++'
    605         write(un_outfor,*) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl)
    606         write(un_outfor,*) dzsnSV(ikl,:)
    607         write(un_outfor,*) TsisSV(ikl,:)
    608         write(un_outfor,*) ro__SV(ikl,:)
    609         write(un_outfor,*) eta_SV(ikl,:)
    610         write(un_outfor,*) G1snSV(ikl,:)
    611         write(un_outfor,*) G2snSV(ikl,:)
    612         write(un_outfor,*) agsnSV(ikl,:)
    613         write(un_outfor,*) istoSV(ikl,:)
    614        
    615       ENDIF
    616 
    617 
    618 
    619 
    620 ! +  -----------------------------                             
    621 ! +  END --- RETURN RESULTS   
    622 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    623       IF (lafin) THEN
    624         fichnom = "restartsis.nc"
    625         CALL sisvatredem("restartsis.nc",ikl2i,rlon,rlat)     
    626        
    627         IF (ok_outfor) THEN
    628           close(unit=un_outfor)                   
    629         END IF
    630       END IF                                                         
    631 
    632 
    633   END SUBROUTINE surf_inlandsis
    634 
    635 
    636 
    637 
    638 
    639 
    640 
    641 
    642 
    643 
    644 
    645 
    646 
    647 
    648 !=======================================================================
    649 
    650   SUBROUTINE get_soil_levels(dz1, dz2, lambda)
    651 ! ======================================================================
    652 ! Routine to compute the vertical discretization of the soil in analogy
    653 ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
    654 ! of SISVAT, therefore it's needed here.
    655 !
    656     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    657     USE mod_phys_lmdz_para
    658     USE VAR_SV
    659 
    660 
    661 !    INCLUDE "dimsoil.h"
    662 
    663     REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
    664     REAL, INTENT(OUT)                     :: lambda
    665 
    666 
    667 !-----------------------------------------------------------------------
    668 !   Depthts:
    669 !   --------
    670     REAL fz,rk,fz1,rk1,rk2
    671     REAL min_period, dalph_soil
    672     INTEGER ierr,jk
    673 
    674     fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
    675 
    676 !    write(*,*)'Start soil level computation'
    677 !-----------------------------------------------------------------------
    678 ! Calculation of some constants
    679 ! NB! These constants do not depend on the sub-surfaces
    680 !-----------------------------------------------------------------------
    681 !-----------------------------------------------------------------------
    682 !   ground levels
    683 !   grnd=z/l where l is the skin depth of the diurnal cycle:
    684 !-----------------------------------------------------------------------
    685 
    686      min_period=1800. ! en secondes
    687      dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
    688 ! !$OMP MASTER
    689 !     IF (is_mpi_root) THEN
    690 !        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
    691 !        IF (ierr == 0) THEN ! Read file only if it exists
    692 !           READ(99,*) min_period
    693 !           READ(99,*) dalph_soil
    694 !           PRINT*,'Discretization for the soil model'
    695 !           PRINT*,'First level e-folding depth',min_period, &
    696 !                '   dalph',dalph_soil
    697 !           CLOSE(99)
    698 !        END IF
    699 !     ENDIF
    700 ! !$OMP END MASTER
    701 !     CALL bcast(min_period)
    702 !     CALL bcast(dalph_soil)
    703 
    704 !   la premiere couche represente un dixieme de cycle diurne
    705      fz1=SQRT(min_period/3.14)
    706      
    707      DO jk=1,nsoilmx
    708         rk1=jk
    709         rk2=jk-1
    710         dz2(jk)=fz(rk1)-fz(rk2)
    711      ENDDO
    712      DO jk=1,nsoilmx-1
    713         rk1=jk+.5
    714         rk2=jk-.5
    715         dz1(jk)=1./(fz(rk1)-fz(rk2))
    716      ENDDO
    717      lambda=fz(.5)*dz1(1)
    718      PRINT*,'full layers, intermediate layers (seconds)'
    719      DO jk=1,nsoilmx
    720         rk=jk
    721         rk1=jk+.5
    722         rk2=jk-.5
    723         PRINT *,'fz=', &
    724              fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
    725      ENDDO
    726 
    727   END SUBROUTINE get_soil_levels
    728  
    729 
    730 
    731 
    732 
    733 
    734 
    735 
    736 
    737 
    738 
    739 
    740 
    741 
    742 
    743 
    744 
    745 !===========================================================================
    746 
    747   SUBROUTINE SISVAT_ini(knon)                                                     
    748                                                                                
    749 !C +------------------------------------------------------------------------+ 
    750 !C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR |
    751 !C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters | 
    752 !C +------------------------------------------------------------------------+ 
    753 !C |   PARAMETERS:  klonv: Total Number of columns =                        | 
    754 !C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       | 
    755 !C |                     X       Number of Mosaic Cell per grid box         | 
    756 !C |                                                                        | 
    757 !C |   INPUT:   dt__SV   : Time  Step                                   [s] |
    758 !C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] | 
    759 !C |                                                                        | 
    760 !C |   OUTPUT:             [-] | 
    761 !C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] | 
    762 !C |            etamSV   : Soil Minimum Humidity                    [m3/m3] | 
    763 !C |                      (based on a prescribed Soil Relative Humidity)    | 
    764 !C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       | 
    765 !C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        | 
    766 !C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   | 
    767 !C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] | 
    768 !C |            dzsnSV(0): Soil first Layer Thickness                   [m] | 
    769 !C |            dzmiSV   : Distance between two contiguous levels       [m] | 
    770 !C |            dz78SV   : 7/8 (Layer Thickness)                        [m] |
    771 !C |            dz34SV   : 3/4 (Layer Thickness)                        [m] |
    772 !C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] |
    773 !C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] |
    774 !C |            dtz_SV   : dt/dz                                      [s/m] |
    775 !C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
    776 !C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      | 
    777 !C |            Explic   : Explicit Parameter = 1.0 - Implic                | 
    778 !C |                                                                        |
    779 !C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
    780 !C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
    781 !C | #          #SH: Hapex-Sahel Values                                     ! 
    782 !C |                                                                        |
    783 !C +------------------------------------------------------------------------+ 
    784 !                                                                             
    785 !                                                                             
    786                                                                              
    787 !C +--Global Variables                                                         
    788 !C +  ================         
    789 
    790       USE dimphy
    791       USE VARphy                                           
    792       USE VAR_SV                                                     
    793       USE VARdSV                                                         
    794       USE VAR0SV                                                           
    795       USE VARxSV
    796       USE VARtSV
    797       USE VARxSV
    798       USE VARySV
    799       IMPLICIT NONE                                                           
    800                                                                              
    801                                                                              
    802                                                                              
    803 !C +--Arguments                                                     
    804 !C +  ==================                                                     
    805        INTEGER,INTENT(IN) ::  knon                                       
    806 
    807 !C +--Internal Variables                                                     
    808 !C +  ==================                                                     
    809                                                                                
    810       INTEGER ::  ivt   ,ist   ,ikl   ,isl   ,isn   ,ikh               
    811       INTEGER ::  misl_2,nisl_2                                             
    812       REAL    ::  d__eta,eta__1,eta__2,Khyd_1,Khyd_2                           
    813       REAL,PARAMETER  ::  RHsMin=  0.001        ! Min.Soil Relative Humidity   
    814       REAL    ::  PsiMax                        ! Max.Soil Water    Potential
    815       REAL    ::  a_Khyd,b_Khyd                 ! Piecewis.https://www.lequipe.fr/Water Conductivity
    816 
    817 
    818 !c #WR REAL    ::  Khyd_x,Khyd_y                                               
    819                                                                              
    820                                                                              
    821                                                                  
    822 !C +--Non Time Dependant SISVAT parameters                                   
    823 !C +  ====================================                               
    824                                                                              
    825 !C +--Soil Discretization                                                     
    826 !C +  -------------------                                                     
    827                                                                              
    828 !C +--Numerical Scheme Parameters                                             
    829 !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^                                             
     742
     743    END SUBROUTINE surf_inlandsis
     744
     745
     746    !=======================================================================
     747
     748    SUBROUTINE get_soil_levels(dz1, dz2, lambda)
     749        ! ======================================================================
     750        ! Routine to compute the vertical discretization of the soil in analogy
     751        ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
     752        ! of SISVAT, therefore it's needed here.
     753        !
     754        USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
     755        USE mod_phys_lmdz_para
     756        USE VAR_SV
     757
     758
     759        !    INCLUDE "dimsoil.h"
     760
     761        REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
     762        REAL, INTENT(OUT) :: lambda
     763
     764
     765        !-----------------------------------------------------------------------
     766        !   Depthts:
     767        !   --------
     768        REAL fz, rk, fz1, rk1, rk2
     769        REAL min_period, dalph_soil
     770        INTEGER ierr, jk
     771
     772        fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
     773
     774        !    write(*,*)'Start soil level computation'
     775        !-----------------------------------------------------------------------
     776        ! Calculation of some constants
     777        ! NB! These constants do not depend on the sub-surfaces
     778        !-----------------------------------------------------------------------
     779        !-----------------------------------------------------------------------
     780        !   ground levels
     781        !   grnd=z/l where l is the skin depth of the diurnal cycle:
     782        !-----------------------------------------------------------------------
     783
     784        min_period = 1800. ! en secondes
     785        dalph_soil = 2.    ! rapport entre les epaisseurs de 2 couches succ.
     786        ! !$OMP MASTER
     787        !     IF (is_mpi_root) THEN
     788        !        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
     789        !        IF (ierr == 0) THEN ! Read file only if it exists
     790        !           READ(99,*) min_period
     791        !           READ(99,*) dalph_soil
     792        !           PRINT*,'Discretization for the soil model'
     793        !           PRINT*,'First level e-folding depth',min_period, &
     794        !                '   dalph',dalph_soil
     795        !           CLOSE(99)
     796        !        END IF
     797        !     ENDIF
     798        ! !$OMP END MASTER
     799        !     CALL bcast(min_period)
     800        !     CALL bcast(dalph_soil)
     801
     802        !   la premiere couche represente un dixieme de cycle diurne
     803        fz1 = SQRT(min_period / 3.14)
     804
     805        DO jk = 1, nsoilmx
     806            rk1 = jk
     807            rk2 = jk - 1
     808            dz2(jk) = fz(rk1) - fz(rk2)
     809        ENDDO
     810        DO jk = 1, nsoilmx - 1
     811            rk1 = jk + .5
     812            rk2 = jk - .5
     813            dz1(jk) = 1. / (fz(rk1) - fz(rk2))
     814        ENDDO
     815        lambda = fz(.5) * dz1(1)
     816        DO jk = 1, nsoilmx
     817            rk = jk
     818            rk1 = jk + .5
     819            rk2 = jk - .5
     820        ENDDO
     821
     822    END SUBROUTINE get_soil_levels
     823
     824
     825    !===========================================================================
     826
     827    SUBROUTINE SISVAT_ini(knon)
     828
     829        !C +------------------------------------------------------------------------+
     830        !C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR |
     831        !C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters |
     832        !C +------------------------------------------------------------------------+
     833        !C |   PARAMETERS:  klonv: Total Number of columns =                        |
     834        !C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
     835        !C |                     X       Number of Mosaic Cell per grid box         |
     836        !C |                                                                        |
     837        !C |   INPUT:   dt__SV   : Time  Step                                   [s] |
     838        !C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] |
     839        !C |                                                                        |
     840        !C |   OUTPUT:             [-] |
     841        !C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] |
     842        !C |            etamSV   : Soil Minimum Humidity                    [m3/m3] |
     843        !C |                      (based on a prescribed Soil Relative Humidity)    |
     844        !C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       |
     845        !C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        |
     846        !C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   |
     847        !C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] |
     848        !C |            dzsnSV(0): Soil first Layer Thickness                   [m] |
     849        !C |            dzmiSV   : Distance between two contiguous levels       [m] |
     850        !C |            dz78SV   : 7/8 (Layer Thickness)                        [m] |
     851        !C |            dz34SV   : 3/4 (Layer Thickness)                        [m] |
     852        !C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] |
     853        !C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] |
     854        !C |            dtz_SV   : dt/dz                                      [s/m] |
     855        !C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
     856        !C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      |
     857        !C |            Explic   : Explicit Parameter = 1.0 - Implic                |
     858        !C |                                                                        |
     859        !C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
     860        !C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
     861        !C | #          #SH: Hapex-Sahel Values                                     !
     862        !C |                                                                        |
     863        !C +------------------------------------------------------------------------+
     864        !
     865        !
     866
     867        !C +--Global Variables
     868        !C +  ================
     869
     870        USE dimphy
     871        USE VARphy
     872        USE VAR_SV
     873        USE VARdSV
     874        USE VAR0SV
     875        USE VARxSV
     876        USE VARtSV
     877        USE VARxSV
     878        USE VARySV
     879        IMPLICIT NONE
     880
     881
     882
     883        !C +--Arguments
     884        !C +  ==================
     885        INTEGER, INTENT(IN) :: knon
     886
     887        !C +--Internal Variables
     888        !C +  ==================
     889
     890        INTEGER :: ivt, ist, ikl, isl, isn, ikh
     891        INTEGER :: misl_2, nisl_2
     892        REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2
     893        REAL, PARAMETER :: RHsMin = 0.001        ! Min.Soil Relative Humidity
     894        REAL :: PsiMax                        ! Max.Soil Water    Potential
     895        REAL :: a_Khyd, b_Khyd                 ! Water conductivity
     896
     897
     898        !c #WR REAL    ::  Khyd_x,Khyd_y
     899
     900
     901
     902        !C +--Non Time Dependant SISVAT parameters
     903        !C +  ====================================
     904
     905        !C +--Soil Discretization
     906        !C +  -------------------
     907
     908        !C +--Numerical Scheme Parameters
     909        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
    830910        Implic = 0.75                           ! 0.5  <==> Crank-Nicholson 
    831911        Explic = 1.00 - Implic                  !                           
    832                                                                              
    833 !C +--Soil/Snow Layers Indices                                               
    834 !C +  ^^^^^^^^^^^^^^^^^^^^^^^^                                               
    835       DO  isl=-nsol,0                                                         
    836         islpSV(isl) =           isl+1                                         
    837         islpSV(isl) = min(      islpSV(isl),0)                               
    838         islmSV(isl) =           isl-1                                         
    839         islmSV(isl) = max(-nsol,islmSV(isl))                                 
    840       END DO                                                                 
    841                                                                                
    842       DO  isn=1,nsno                                                           
    843         isnpSV(isn) =           isn+1                                         
    844         isnpSV(isn) = min(      isnpSV(isn),nsno)                           
    845       END DO                                                                 
    846                                                                              
    847 !C +--Soil      Layers Thicknesses                                             
    848 !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
    849 ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!   
    850 !c #kd IF (nsol.gt.4)                                              THEN       
    851 !c #kd   DO isl=-5,-nsol,-1                                                   
    852 !c #kd     dz_dSV(isl)=   1.                                                 
    853 !c #kd   END DO                                                               
    854 !c #kd END IF                                                                 
    855 !                                                                             
    856 !      IF (nsol.ne.4)                                              THEN       
    857 !        DO isl= 0,-nsol,-1                                                   
    858 !          misl_2 =     -mod(isl,2)                                         
    859 !          nisl_2 =         -isl/2                                           
    860 !          dz_dSV(isl)=(((1-misl_2) * 0.001                                   
    861 !     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.             
    862 !C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm           
    863 !                                                                             
    864 !c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001                                   
    865 !c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.             
    866 !                                                                             
    867 !c #05     dz_dSV(isl)=(((1-misl_2) * 0.001                                   
    868 !c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5             
    869 !        END DO                                                               
    870 !          dz_dSV(0)  =               0.001                                   
    871 !          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004         
    872 !      END IF           
    873 
    874                                                                                
    875         zz_dSV      = 0.                                                     
    876       DO  isl=-nsol,0                                                       
    877         dzmiSV(isl) = 0.500*(dz_dSV(isl)        +dz_dSV(islmSV(isl)))       
    878         dziiSV(isl) = 0.500* dz_dSV(isl)        /dzmiSV(isl)                 
    879         dzi_SV(isl) = 0.500* dz_dSV(islmSV(isl))/dzmiSV(isl)                   
    880         dtz_SV(isl) =        dt__SV             /dz_dSV(isl) 
    881         dtz_SV2(isl) =        1.            /dz_dSV(isl)                         
    882         dz78SV(isl) = 0.875* dz_dSV(isl)                                     
    883         dz34SV(isl) = 0.750* dz_dSV(isl)                                     
    884         dz_8SV(isl) = 0.125* dz_dSV(isl)                                       
    885         dzAvSV(isl) = 0.125* dz_dSV(islmSV(isl))                        &
    886      &              + 0.750* dz_dSV(isl)                                &
    887      &              + 0.125* dz_dSV(islpSV(isl))                                                             
    888         zz_dSV      = zz_dSV+dz_dSV(isl)                                     
    889       END DO                                                                 
    890       DO ikl=1,knon !v                                                         
    891         dzsnSV(ikl,0) =      dz_dSV(0)                                       
    892       END DO                                                                 
    893                                                                              
    894 !C +--Conversion to a 50 m Swab Ocean Discretization                           
    895 !C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                           
    896         OcndSV = 0.                                                           
    897       DO isl=-nsol,0                                                         
    898         OcndSV = OcndSV +dz_dSV(isl)                                           
    899       END DO                                                                   
    900         OcndSV = 50.    /OcndSV                                               
    901                                                                                                            
    902                                                                                                                                                                                                                                
    903 !C +--Secondary Soil       Parameters                                         
    904 !C +  -------------------------------                                         
    905                                                                                
    906       DO  ist=0,nsot                                                           
    907          rocsSV(ist)=(1.0-etadSV(ist))*1.2E+6   ! Soil Contrib. to (ro c)_s   
    908          s1__SV(ist)=     bCHdSV(ist)          & ! Factor of (eta)**(b+2)     
    909      &  *psidSV(ist)     *Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)     
    910      & /(etadSV(ist)**(   bCHdSV(ist)+3.))     !                             
    911          s2__SV(ist)=     Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)     
    912      & /(etadSV(ist)**(2.*bCHdSV(ist)+3.))     !    in DR97, Eqn.(3.35)       
    913                                                                              
    914 !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)     
    915 !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^   
    916          Psimax = -(log(RHsMin))/7.2E-5        ! DR97, Eqn 3.15 Inversion   
    917          etamSV(ist) =  etadSV(ist)                                      &
    918      &         *(PsiMax/psidSV(ist))**(-min(10.,1./bCHdSV(ist)))             
    919       END DO                                                                 
    920          etamSV(12)  =  0.                                                   
    921                                                                              
    922 !C +--Piecewise Hydraulic Conductivity Profiles                               
    923 !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
    924       DO   ist=0,nsot                                                         
    925  
    926                                                                              
    927           d__eta          =  etadSV(ist)/nkhy                                 
    928           eta__1          =  0.                                             
    929           eta__2          =  d__eta                                           
    930         DO ikh=0,nkhy                                                         
    931           Khyd_1          =  s2__SV(ist)             & ! DR97, Eqn.(3.35)     
    932      &  *(eta__1      **(2. *bCHdSV(ist)+3.))        !                       
    933           Khyd_2          =  s2__SV(ist)             &!                       
    934      &  *(eta__2      **(2. *bCHdSV(ist)+3.))        !                       
    935                                                                              
    936           a_Khyd          = (Khyd_2-Khyd_1)/d__eta   !                       
    937           b_Khyd          =  Khyd_1-a_Khyd *eta__1   !                       
    938 !c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !                       
    939 !c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !                       
    940           aKdtSV(ist,ikh) =  a_Khyd       * dt__SV   !                       
    941           bKdtSV(ist,ikh) =  b_Khyd       * dt__SV   !                         
    942                
    943           eta__1          = eta__1  + d__eta                             
    944           eta__2          = eta__2  + d__eta                             
    945         END DO                                                             
    946       END DO                                                               
    947 
    948                                                                            
    949       return                                                               
    950 
    951   END SUBROUTINE SISVAT_ini
    952 
    953 
    954 
    955 
    956 
    957 
    958 
    959 !***************************************************************************
    960 
    961     SUBROUTINE sisvatetat0 (fichnom,ikl2i)
    962 
    963     USE dimphy
    964     USE mod_grid_phy_lmdz
    965     USE mod_phys_lmdz_para
    966 
    967     USE iostart
    968     USE VAR_SV
    969     USE VARdSV
    970     USE VARxSV         
    971     USE VARtSV
    972     USE indice_sol_mod
    973 
    974       IMPLICIT none
    975 !======================================================================
    976 ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
    977 ! Objet: Lecture du fichier de conditions initiales pour SISVAT
    978 !======================================================================
    979     include "netcdf.inc"
    980 !    include "indicesol.h"
    981 
    982 !    include "dimsoil.h"
    983     include "clesphys.h"
    984     include "thermcell.h"
    985     include "compbl.h"
    986 
    987 !======================================================================
    988     CHARACTER(LEN=*) :: fichnom
    989 
    990 
    991     INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
    992     REAL, DIMENSION(klon) :: rlon
    993     REAL, DIMENSION(klon) :: rlat
    994 
    995 ! les variables globales ecrites dans le fichier restart
    996     REAL, DIMENSION(klon) :: isno
    997     REAL, DIMENSION(klon) :: ispi
    998     REAL, DIMENSION(klon) :: iice
    999     REAL, DIMENSION(klon) :: rusn
    1000     REAL, DIMENSION(klon, nsno) :: isto
    1001 
    1002     REAL, DIMENSION(klon, nsismx) :: Tsis
    1003     REAL, DIMENSION(klon, nsismx) :: eta
    1004     REAL, DIMENSION(klon, nsismx) :: ro
    1005 
    1006     REAL, DIMENSION(klon, nsno) :: dzsn     
    1007     REAL, DIMENSION(klon, nsno) :: G1sn
    1008     REAL, DIMENSION(klon, nsno) :: G2sn
    1009     REAL, DIMENSION(klon, nsno) :: agsn
    1010 
    1011     REAL, DIMENSION(klon) :: toic
    1012 
    1013 
    1014     INTEGER  :: isl, ikl, i, isn , errT, erreta, errro, errdz, snopts
    1015     CHARACTER (len=2) :: str2
    1016     LOGICAL :: found
    1017  
    1018     errT=0
    1019     errro=0
    1020     erreta=0
    1021     errdz=0
    1022     snopts=0
    1023 ! Ouvrir le fichier contenant l'etat initial:
    1024 
    1025       CALL open_startphy(fichnom)
    1026 
    1027 ! Lecture des latitudes, longitudes (coordonnees):
    1028 
    1029       CALL get_field("latitude",rlat,found)
    1030       CALL get_field("longitude",rlon,found)
    1031 
    1032       CALL get_field("n_snows", isno,found)
    1033       IF (.NOT. found) THEN
    1034         PRINT*, 'phyetat0: Le champ <n_snows> est absent'
    1035         PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
    1036       ENDIF
    1037 
    1038       CALL get_field("n_ice_top",ispi,found)
    1039       CALL get_field("n_ice",iice,found)
    1040       CALL get_field("surf_water",rusn,found)
    1041 !      IF (.NOT. found) THEN
    1042 !        PRINT*, 'phyetat0: Le champ <surf_water> est absent'
    1043 !        rusn(:)=0. 
    1044 !      ENDIF
    1045 
    1046 
    1047       CALL get_field("to_ice",toic,found)
    1048       IF (.NOT. found) THEN
    1049         PRINT*, 'phyetat0: Le champ <to_ice> est absent'
    1050         toic(:)=0. 
    1051       ENDIF
    1052 
    1053 
    1054 
    1055       DO isn = 1,nsno
    1056         IF (isn.LE.99) THEN
    1057             WRITE(str2,'(i2.2)') isn
    1058             CALL get_field("AGESNOW"//str2, &
    1059                           agsn(:,isn),found)
    1060         ELSE
    1061             PRINT*, "Trop de couches"
    1062             CALL abort
     912
     913        !C +--Soil/Snow Layers Indices
     914        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^
     915        DO  isl = -nsol, 0
     916            islpSV(isl) = isl + 1
     917            islpSV(isl) = min(islpSV(isl), 0)
     918            islmSV(isl) = isl - 1
     919            islmSV(isl) = max(-nsol, islmSV(isl))
     920        END DO
     921
     922        DO  isn = 1, nsno
     923            isnpSV(isn) = isn + 1
     924            isnpSV(isn) = min(isnpSV(isn), nsno)
     925        END DO
     926
     927        !C +--Soil      Layers Thicknesses
     928        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     929        ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
     930        !c #kd IF (nsol.gt.4)                                              THEN
     931        !c #kd   DO isl=-5,-nsol,-1
     932        !c #kd     dz_dSV(isl)=   1.
     933        !c #kd   END DO
     934        !c #kd END IF
     935        !
     936        !      IF (nsol.ne.4)                                              THEN
     937        !        DO isl= 0,-nsol,-1
     938        !          misl_2 =     -mod(isl,2)
     939        !          nisl_2 =         -isl/2
     940        !          dz_dSV(isl)=(((1-misl_2) * 0.001
     941        !     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.
     942        !C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm
     943        !
     944        !c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001
     945        !c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.
     946        !
     947        !c #05     dz_dSV(isl)=(((1-misl_2) * 0.001
     948        !c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5
     949        !        END DO
     950        !          dz_dSV(0)  =               0.001
     951        !          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004
     952        !      END IF
     953
     954        zz_dSV = 0.
     955        DO  isl = -nsol, 0
     956            dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl)))
     957            dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl)
     958            dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl)
     959            dtz_SV(isl) = dt__SV / dz_dSV(isl)
     960            dtz_SV2(isl) = 1. / dz_dSV(isl)
     961            dz78SV(isl) = 0.875 * dz_dSV(isl)
     962            dz34SV(isl) = 0.750 * dz_dSV(isl)
     963            dz_8SV(isl) = 0.125 * dz_dSV(isl)
     964            dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl))                        &
     965                    & + 0.750 * dz_dSV(isl)                                &
     966                    & + 0.125 * dz_dSV(islpSV(isl))
     967            zz_dSV = zz_dSV + dz_dSV(isl)
     968        END DO
     969        DO ikl = 1, knon !v
     970            dzsnSV(ikl, 0) = dz_dSV(0)
     971        END DO
     972
     973        !C +--Conversion to a 50 m Swab Ocean Discretization
     974        !C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     975        OcndSV = 0.
     976        DO isl = -nsol, 0
     977            OcndSV = OcndSV + dz_dSV(isl)
     978        END DO
     979        OcndSV = 50. / OcndSV
     980
     981
     982        !C +--Secondary Soil       Parameters
     983        !C +  -------------------------------
     984
     985        DO  ist = 0, nsot
     986            rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6   ! Soil Contrib. to (ro c)_s
     987            s1__SV(ist) = bCHdSV(ist)          & ! Factor of (eta)**(b+2)
     988                    & * psidSV(ist) * Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)
     989                    & / (etadSV(ist)**(bCHdSV(ist) + 3.))     !
     990            s2__SV(ist) = Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)
     991                    & / (etadSV(ist)**(2. * bCHdSV(ist) + 3.))     !    in DR97, Eqn.(3.35)
     992
     993            !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
     994            !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     995            Psimax = -(log(RHsMin)) / 7.2E-5        ! DR97, Eqn 3.15 Inversion
     996            etamSV(ist) = etadSV(ist)                                      &
     997                    & * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist)))
     998        END DO
     999        etamSV(12) = 0.
     1000
     1001        !C +--Piecewise Hydraulic Conductivity Profiles
     1002        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     1003        DO   ist = 0, nsot
     1004
     1005            d__eta = etadSV(ist) / nkhy
     1006            eta__1 = 0.
     1007            eta__2 = d__eta
     1008            DO ikh = 0, nkhy
     1009                Khyd_1 = s2__SV(ist)             & ! DR97, Eqn.(3.35)
     1010                        & * (eta__1      **(2. * bCHdSV(ist) + 3.))        !
     1011                Khyd_2 = s2__SV(ist)             &!
     1012                        & * (eta__2      **(2. * bCHdSV(ist) + 3.))        !
     1013
     1014                a_Khyd = (Khyd_2 - Khyd_1) / d__eta   !
     1015                b_Khyd = Khyd_1 - a_Khyd * eta__1   !
     1016                !c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !
     1017                !c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !
     1018                aKdtSV(ist, ikh) = a_Khyd * dt__SV   !
     1019                bKdtSV(ist, ikh) = b_Khyd * dt__SV   !
     1020
     1021                eta__1 = eta__1 + d__eta
     1022                eta__2 = eta__2 + d__eta
     1023            END DO
     1024        END DO
     1025
     1026        return
     1027
     1028    END SUBROUTINE SISVAT_ini
     1029
     1030
     1031    !***************************************************************************
     1032
     1033    SUBROUTINE sisvatetat0 (fichnom, ikl2i)
     1034
     1035        USE dimphy
     1036        USE mod_grid_phy_lmdz
     1037        USE mod_phys_lmdz_para
     1038
     1039        USE iostart
     1040        USE VAR_SV
     1041        USE VARdSV
     1042        USE VARxSV
     1043        USE VARtSV
     1044        USE indice_sol_mod
     1045
     1046        IMPLICIT none
     1047        !======================================================================
     1048        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
     1049        ! Objet: Lecture du fichier de conditions initiales pour SISVAT
     1050        !======================================================================
     1051        include "netcdf.inc"
     1052        !    include "indicesol.h"
     1053
     1054        !    include "dimsoil.h"
     1055        include "clesphys.h"
     1056        include "thermcell.h"
     1057        include "compbl.h"
     1058
     1059        !======================================================================
     1060        CHARACTER(LEN = *) :: fichnom
     1061
     1062        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
     1063        REAL, DIMENSION(klon) :: rlon
     1064        REAL, DIMENSION(klon) :: rlat
     1065
     1066        ! les variables globales ecrites dans le fichier restart
     1067        REAL, DIMENSION(klon) :: isno
     1068        REAL, DIMENSION(klon) :: ispi
     1069        REAL, DIMENSION(klon) :: iice
     1070        REAL, DIMENSION(klon) :: rusn
     1071        REAL, DIMENSION(klon, nsno) :: isto
     1072
     1073        REAL, DIMENSION(klon, nsismx) :: Tsis
     1074        REAL, DIMENSION(klon, nsismx) :: eta
     1075        REAL, DIMENSION(klon, nsismx) :: ro
     1076
     1077        REAL, DIMENSION(klon, nsno) :: dzsn
     1078        REAL, DIMENSION(klon, nsno) :: G1sn
     1079        REAL, DIMENSION(klon, nsno) :: G2sn
     1080        REAL, DIMENSION(klon, nsno) :: agsn
     1081
     1082        REAL, DIMENSION(klon) :: toic
     1083
     1084        INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts
     1085        CHARACTER (len = 2) :: str2
     1086        LOGICAL :: found
     1087
     1088        errT = 0
     1089        errro = 0
     1090        erreta = 0
     1091        errdz = 0
     1092        snopts = 0
     1093        ! Ouvrir le fichier contenant l'etat initial:
     1094
     1095        CALL open_startphy(fichnom)
     1096
     1097        ! Lecture des latitudes, longitudes (coordonnees):
     1098
     1099        CALL get_field("latitude", rlat, found)
     1100        CALL get_field("longitude", rlon, found)
     1101
     1102        CALL get_field("n_snows", isno, found)
     1103        IF (.NOT. found) THEN
     1104            PRINT*, 'phyetat0: Le champ <n_snows> est absent'
     1105            PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
    10631106        ENDIF
    1064       ENDDO
    1065       DO isn = 1,nsno
    1066         IF (isn.LE.99) THEN
    1067             WRITE(str2,'(i2.2)') isn
    1068             CALL get_field("DZSNOW"//str2, &
    1069                           dzsn(:,isn),found)
    1070         ELSE
    1071             PRINT*, "Trop de couches"
    1072             CALL abort
     1107
     1108        CALL get_field("n_ice_top", ispi, found)
     1109        CALL get_field("n_ice", iice, found)
     1110        CALL get_field("surf_water", rusn, found)
     1111
     1112
     1113        CALL get_field("to_ice", toic, found)
     1114        IF (.NOT. found) THEN
     1115            PRINT*, 'phyetat0: Le champ <to_ice> est absent'
     1116            toic(:) = 0.
    10731117        ENDIF
    1074       ENDDO
    1075       DO isn = 1,nsno
    1076         IF (isn.LE.99) THEN
    1077             WRITE(str2,'(i2.2)') isn
    1078             CALL get_field("G2SNOW"//str2, &
    1079                           G2sn(:,isn),found)
    1080         ELSE
    1081             PRINT*, "Trop de couches"
    1082             CALL abort
    1083         ENDIF
    1084       ENDDO
    1085       DO isn = 1,nsno
    1086         IF (isn.LE.99) THEN
    1087             WRITE(str2,'(i2.2)') isn
    1088             CALL get_field("G1SNOW"//str2, &
    1089                           G1sn(:,isn),found)
    1090         ELSE
    1091             PRINT*, "Trop de couches"
    1092             CALL abort
    1093         ENDIF
    1094       ENDDO
    1095       DO isn = 1,nsismx
    1096         IF (isn.LE.99) THEN
    1097             WRITE(str2,'(i2.2)') isn
    1098             CALL get_field("ETA"//str2, &
    1099                           eta(:,isn),found)
    1100         ELSE
    1101             PRINT*, "Trop de couches"
    1102             CALL abort
    1103         ENDIF
    1104       ENDDO
    1105       DO isn = 1,nsismx
    1106         IF (isn.LE.99) THEN
    1107             WRITE(str2,'(i2.2)') isn
    1108             CALL get_field("RO"//str2, &
    1109                           ro(:,isn),found)
    1110         ELSE
    1111             PRINT*, "Trop de couches"
    1112             CALL abort
    1113         ENDIF
    1114       ENDDO
    1115       DO isn = 1,nsismx
    1116         IF (isn.LE.99) THEN
    1117             WRITE(str2,'(i2.2)') isn
    1118             CALL get_field("TSS"//str2, &
    1119                           Tsis(:,isn),found)
    1120         ELSE
    1121             PRINT*, "Trop de couches"
    1122             CALL abort
    1123         ENDIF
    1124       ENDDO
    1125       DO isn = 1,nsno
    1126         IF (isn.LE.99) THEN
    1127             WRITE(str2,'(i2.2)') isn
    1128             CALL get_field("HISTORY"//str2, &
    1129                           isto(:,isn),found)
    1130         ELSE
    1131             PRINT*, "Trop de couches"
    1132             CALL abort
    1133         ENDIF
    1134       ENDDO
    1135       write(*,*)'Read ',fichnom,' finished!!'
    1136 
    1137 !*********************************************************************************
    1138 ! Compress restart file variables for SISVAT
    1139 
    1140 
    1141       DO  ikl = 1,klon                                                   
    1142           i   = ikl2i(ikl)   
    1143           IF (i > 0) THEN
    1144               isnoSV(ikl)     = INT(isno(i))          ! Nb Snow/Ice Lay.   
    1145               ispiSV(ikl)     = INT(ispi(i))          ! Nb Supr.Ice Lay. 
    1146               iiceSV(ikl)     = INT(iice(i))          ! Nb      Ice Lay.   
    1147                                                                              
    1148 
    1149             DO isl =   -nsol,0   
    1150               ro__SV(ikl,isl) = ro(i,nsno+1-isl)       !                   
    1151               eta_SV(ikl,isl) = eta(i,nsno+1-isl)         ! Soil Humidity     
    1152 !hjp 15/10/2010
    1153               IF (eta_SV(ikl,isl) <= 1.e-6) THEN          !hj check
    1154                 eta_SV(ikl,isl) = 1.e-6
    1155               ENDIF
    1156               TsisSV(ikl,isl) = Tsis(i,nsno+1-isl)        ! Soil Temperature 
    1157               IF (TsisSV(ikl,isl) <= 1.) THEN             !hj check
    1158 !                errT=errT+1
    1159                 TsisSV(ikl,isl) = 273.15-0.2              ! Etienne: negative temperature since soil is ice
    1160               ENDIF
    1161 
    1162             END DO       
    1163             write(*,*)'Copy histo', ikl
    1164    
    1165    
    1166             DO  isn = 1,isnoSV(ikl) !nsno     
    1167               snopts=snopts+1
    1168               IF (isto(i,isn) > 10.) THEN          !hj check
    1169                 write(*,*)'Irregular isto',ikl,i,isn,isto(i,isn)
    1170                 isto(i,isn) = 1.
    1171               ENDIF
    1172 
    1173               istoSV(ikl,isn) = INT(isto(i,isn))     ! Snow     History
    1174               ro__SV(ikl,isn) = ro(i,isn)            !        [kg/m3]     
    1175               eta_SV(ikl,isn) = eta(i,isn)           !        [m3/m3] 
    1176               TsisSV(ikl,isn) = Tsis(i,isn)          !            [K] 
    1177 
    1178              IF (TsisSV(ikl,isn) <= 1.) THEN          !hj check
    1179               errT=errT+1
    1180               TsisSV(ikl,isn) = TsisSV(ikl,0)
    1181              ENDIF 
    1182              IF (TsisSV(ikl,isn) <= 1.) THEN          !hj check
    1183               TsisSV(ikl,isn) = 263.15
    1184              ENDIF
    1185              IF (eta_SV(ikl,isn) < 1.e-9) THEN          !hj check
    1186               eta_SV(ikl,isn) = 1.e-6 
    1187               erreta=erreta+1
    1188              ENDIF 
    1189              IF (ro__SV(ikl,isn) <= 10.) THEN          !hj check
    1190               ro__SV(ikl,isn) = 11.
    1191               errro=errro+1
    1192              ENDIF
    1193               write(*,*)ikl,i,isn,Tsis(i,isn),G1sn(i,isn)
    1194               G1snSV(ikl,isn) = G1sn(i,isn)          ! [-]        [-]     
    1195               G2snSV(ikl,isn) = G2sn(i,isn)          ! [-] [0.0001 m]
    1196               dzsnSV(ikl,isn) = dzsn(i,isn)          !            [m]         
    1197               agsnSV(ikl,isn) = agsn(i,isn)          !          [day]     
    1198             END DO 
    1199               rusnSV(ikl)     = rusn(i)              ! Surficial Water   
    1200               toicSV(ikl)     = toic(i)              ! bilan snow to ice   
    1201           END IF                       
    1202         END DO   
     1118
     1119        DO isn = 1, nsno
     1120            IF (isn.LE.99) THEN
     1121                WRITE(str2, '(i2.2)') isn
     1122                CALL get_field("AGESNOW" // str2, &
     1123                        agsn(:, isn), found)
     1124            ELSE
     1125                PRINT*, "Trop de couches"
     1126                CALL abort
     1127            ENDIF
     1128        ENDDO
     1129        DO isn = 1, nsno
     1130            IF (isn.LE.99) THEN
     1131                WRITE(str2, '(i2.2)') isn
     1132                CALL get_field("DZSNOW" // str2, &
     1133                        dzsn(:, isn), found)
     1134            ELSE
     1135                PRINT*, "Trop de couches"
     1136                CALL abort
     1137            ENDIF
     1138        ENDDO
     1139        DO isn = 1, nsno
     1140            IF (isn.LE.99) THEN
     1141                WRITE(str2, '(i2.2)') isn
     1142                CALL get_field("G2SNOW" // str2, &
     1143                        G2sn(:, isn), found)
     1144            ELSE
     1145                PRINT*, "Trop de couches"
     1146                CALL abort
     1147            ENDIF
     1148        ENDDO
     1149        DO isn = 1, nsno
     1150            IF (isn.LE.99) THEN
     1151                WRITE(str2, '(i2.2)') isn
     1152                CALL get_field("G1SNOW" // str2, &
     1153                        G1sn(:, isn), found)
     1154            ELSE
     1155                PRINT*, "Trop de couches"
     1156                CALL abort
     1157            ENDIF
     1158        ENDDO
     1159        DO isn = 1, nsismx
     1160            IF (isn.LE.99) THEN
     1161                WRITE(str2, '(i2.2)') isn
     1162                CALL get_field("ETA" // str2, &
     1163                        eta(:, isn), found)
     1164            ELSE
     1165                PRINT*, "Trop de couches"
     1166                CALL abort
     1167            ENDIF
     1168        ENDDO
     1169        DO isn = 1, nsismx
     1170            IF (isn.LE.99) THEN
     1171                WRITE(str2, '(i2.2)') isn
     1172                CALL get_field("RO" // str2, &
     1173                        ro(:, isn), found)
     1174            ELSE
     1175                PRINT*, "Trop de couches"
     1176                CALL abort
     1177            ENDIF
     1178        ENDDO
     1179        DO isn = 1, nsismx
     1180            IF (isn.LE.99) THEN
     1181                WRITE(str2, '(i2.2)') isn
     1182                CALL get_field("TSS" // str2, &
     1183                        Tsis(:, isn), found)
     1184            ELSE
     1185                PRINT*, "Trop de couches"
     1186                CALL abort
     1187            ENDIF
     1188        ENDDO
     1189        DO isn = 1, nsno
     1190            IF (isn.LE.99) THEN
     1191                WRITE(str2, '(i2.2)') isn
     1192                CALL get_field("HISTORY" // str2, &
     1193                        isto(:, isn), found)
     1194            ELSE
     1195                PRINT*, "Trop de couches"
     1196                CALL abort
     1197            ENDIF
     1198        ENDDO
     1199        write(*, *)'Read ', fichnom, ' finished!!'
     1200
     1201        !*********************************************************************************
     1202        ! Compress restart file variables for SISVAT
     1203
     1204        DO  ikl = 1, klon
     1205            i = ikl2i(ikl)
     1206            IF (i > 0) THEN
     1207                isnoSV(ikl) = INT(isno(i))          ! Nb Snow/Ice Lay.
     1208                ispiSV(ikl) = INT(ispi(i))          ! Nb Supr.Ice Lay.
     1209                iiceSV(ikl) = INT(iice(i))          ! Nb      Ice Lay.
     1210
     1211                DO isl = -nsol, 0
     1212                    ro__SV(ikl, isl) = ro(i, nsno + 1 - isl)       !
     1213                    eta_SV(ikl, isl) = eta(i, nsno + 1 - isl)         ! Soil Humidity
     1214                    !hjp 15/10/2010
     1215                    IF (eta_SV(ikl, isl) <= 1.e-6) THEN          !hj check
     1216                        eta_SV(ikl, isl) = 1.e-6
     1217                    ENDIF
     1218                    TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl)        ! Soil Temperature
     1219                    IF (TsisSV(ikl, isl) <= 1.) THEN             !hj check
     1220                        !                errT=errT+1
     1221                        TsisSV(ikl, isl) = 273.15 - 0.2              ! Etienne: negative temperature since soil is ice
     1222                    ENDIF
     1223
     1224                END DO
     1225                write(*, *)'Copy histo', ikl
     1226
     1227                DO  isn = 1, isnoSV(ikl) !nsno
     1228                    snopts = snopts + 1
     1229                    IF (isto(i, isn) > 10.) THEN          !hj check
     1230                        write(*, *)'Irregular isto', ikl, i, isn, isto(i, isn)
     1231                        isto(i, isn) = 1.
     1232                    ENDIF
     1233
     1234                    istoSV(ikl, isn) = INT(isto(i, isn))     ! Snow     History
     1235                    ro__SV(ikl, isn) = ro(i, isn)            !        [kg/m3]
     1236                    eta_SV(ikl, isn) = eta(i, isn)           !        [m3/m3]
     1237                    TsisSV(ikl, isn) = Tsis(i, isn)          !            [K]
     1238
     1239                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
     1240                        errT = errT + 1
     1241                        TsisSV(ikl, isn) = TsisSV(ikl, 0)
     1242                    ENDIF
     1243                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
     1244                        TsisSV(ikl, isn) = 263.15
     1245                    ENDIF
     1246                    IF (eta_SV(ikl, isn) < 1.e-9) THEN          !hj check
     1247                        eta_SV(ikl, isn) = 1.e-6
     1248                        erreta = erreta + 1
     1249                    ENDIF
     1250                    IF (ro__SV(ikl, isn) <= 10.) THEN          !hj check
     1251                        ro__SV(ikl, isn) = 11.
     1252                        errro = errro + 1
     1253                    ENDIF
     1254                    write(*, *)ikl, i, isn, Tsis(i, isn), G1sn(i, isn)
     1255                    G1snSV(ikl, isn) = G1sn(i, isn)          ! [-]        [-]
     1256                    G2snSV(ikl, isn) = G2sn(i, isn)          ! [-] [0.0001 m]
     1257                    dzsnSV(ikl, isn) = dzsn(i, isn)          !            [m]
     1258                    agsnSV(ikl, isn) = agsn(i, isn)          !          [day]
     1259                END DO
     1260                rusnSV(ikl) = rusn(i)              ! Surficial Water
     1261                toicSV(ikl) = toic(i)              ! bilan snow to ice
     1262            END IF
     1263        END DO
    12031264
    12041265    END SUBROUTINE sisvatetat0
    12051266
    12061267
    1207 
    1208 
    1209 !======================================================================
    1210     SUBROUTINE sisvatredem (fichnom,ikl2i,rlon,rlat)
    1211    
    1212      
    1213    
    1214 !======================================================================
    1215 ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
    1216 ! Objet: Ecriture de l'etat de redemarrage pour SISVAT
    1217 !======================================================================
    1218     USE mod_grid_phy_lmdz
    1219     USE mod_phys_lmdz_para
    1220     USE iostart
    1221     USE VAR_SV
    1222     USE VARxSV         
    1223     USE VARySV !hj tmp 12 03 2010
    1224     USE VARtSV
    1225     USE indice_sol_mod
    1226     USE dimphy
    1227 
    1228 
    1229     IMPLICIT none
    1230 
    1231     include "netcdf.inc"
    1232 !    include "indicesol.h"
    1233 !    include "dimsoil.h"
    1234     include "clesphys.h"
    1235     include "thermcell.h"
    1236     include "compbl.h"
    1237 
    1238 !======================================================================
    1239 
    1240     CHARACTER(LEN=*) :: fichnom
    1241     INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
    1242     REAL, DIMENSION(klon), INTENT(IN) :: rlon
    1243     REAL, DIMENSION(klon), INTENT(IN) :: rlat
    1244 
    1245 ! les variables globales ecrites dans le fichier restart
    1246     REAL, DIMENSION(klon) :: isno
    1247     REAL, DIMENSION(klon) :: ispi
    1248     REAL, DIMENSION(klon) :: iice
    1249     REAL, DIMENSION(klon, nsnowmx) :: isto
    1250 
    1251     REAL, DIMENSION(klon, nsismx) :: Tsis
    1252     REAL, DIMENSION(klon, nsismx) :: eta
    1253     REAL, DIMENSION(klon, nsnowmx) :: dzsn
    1254     REAL, DIMENSION(klon, nsismx) :: ro       
    1255     REAL, DIMENSION(klon, nsnowmx) :: G1sn
    1256     REAL, DIMENSION(klon, nsnowmx) :: G2sn
    1257     REAL, DIMENSION(klon, nsnowmx) :: agsn
    1258     REAL, DIMENSION(klon) :: IRs
    1259     REAL, DIMENSION(klon) :: LMO
    1260     REAL, DIMENSION(klon) :: rusn
    1261     REAL, DIMENSION(klon) :: toic
    1262     REAL, DIMENSION(klon) :: Bufs
    1263     REAL, DIMENSION(klon) :: alb1,alb2,alb3
    1264 
    1265     INTEGER isl, ikl, i, isn, ierr
    1266     CHARACTER (len=2) :: str2
    1267     INTEGER           :: pass
    1268 
    1269       isno(:)       = 0       
    1270       ispi(:)       = 0
    1271       iice(:)       = 0                               
    1272       IRs(:)        = 0.
    1273       LMO(:)        = 0.                             
    1274       eta(:,:)      = 0.     
    1275       Tsis(:,:)     = 0.         
    1276       isto(:,:)     = 0
    1277       ro(:,:)       = 0.       
    1278       G1sn(:,:)     = 0.   
    1279       G2sn(:,:)     = 0.
    1280       dzsn(:,:)     = 0.       
    1281       agsn(:,:)     = 0.
    1282       rusn(:)       = 0.   
    1283       toic(:)       = 0.   
    1284       Bufs(:)       = 0.   
    1285       alb1(:)       = 0.
    1286       alb2(:)       = 0.
    1287       alb3(:)       = 0.
    1288 
    1289 !***************************************************************************
    1290 ! Uncompress SISVAT output variables for storage
    1291            
    1292 
    1293       print*, 'je rentre dans restart inlandsis'     
    1294       DO  ikl = 1,klon 
    1295            i   = ikl2i(ikl)
    1296       IF (i > 0) THEN
    1297         isno(i)       = 1.*isnoSV(ikl)               ! Nb Snow/Ice Lay.   
    1298         ispi(i)       = 1.*ispiSV(ikl)               ! Nb Supr.Ice Lay.   
    1299         iice(i)       = 1.*iiceSV(ikl)               ! Nb      Ice Lay.       
    1300    
    1301 !        IRs(i)        = IRs_SV(ikl)
    1302 !        LMO(i)        = LMO_SV(ikl)                             
    1303 
    1304 
    1305         DO isl =   -nsol,0                           !                   
    1306           eta(i,nsno+1-isl)  = eta_SV(ikl,isl)            ! Soil Humidity     
    1307           Tsis(i,nsno+1-isl) = TsisSV(ikl,isl)            ! Soil Temperature   
    1308           ro(i,nsno+1-isl)   = ro__SV(ikl,isl)            !        [kg/m3]   
    1309         END DO       
    1310  
    1311  
    1312         DO  isn = 1,nsno             
    1313           isto(i,isn)   = 1.*istoSV(ikl,isn)         ! Snow     History
    1314           ro(i,isn)     = ro__SV(ikl,isn)            !        [kg/m3]     
    1315           eta(i,isn)    = eta_SV(ikl,isn)            !        [m3/m3] 
    1316           Tsis(i,isn)   = TsisSV(ikl,isn)            !            [K]   
    1317           G1sn(i,isn)   = G1snSV(ikl,isn)            ! [-]        [-]     
    1318           G2sn(i,isn)   = G2snSV(ikl,isn)            ! [-] [0.0001 m]
    1319           dzsn(i,isn)   = dzsnSV(ikl,isn)            !            [m]         
    1320           agsn(i,isn)   = agsnSV(ikl,isn)            !          [day]     
    1321         END DO 
    1322         rusn(i)       = rusnSV(ikl)                  ! Surficial Water 
    1323         toic(i)       = toicSV(ikl)                  ! to ice
    1324         alb1(i)       = alb1sv(ikl)
    1325         alb2(i)       = alb2sv(ikl)
    1326         alb3(i)       = alb3sv(ikl)
    1327 !        Bufs(i)       = BufsSV(ikl)     
    1328       END IF                     
    1329       END DO                                               
    1330 
    1331 
    1332       print*, 'je call open_restart'     
    1333 
    1334       CALL open_restartphy(fichnom)
    1335 
    1336       print*, 'je sors open_restart'     
    1337 
    1338 
    1339       DO pass = 1, 2
    1340         CALL put_field(pass,"longitude", &
    1341                     "Longitudes de la grille physique",rlon)     
    1342         CALL put_field(pass,"latitude","Latitudes de la grille physique",rlat)
    1343  
    1344         CALL put_field(pass,"n_snows", "number of snow/ice layers",isno)
    1345         CALL put_field(pass,"n_ice_top", "number of top ice layers",ispi)
    1346         CALL put_field(pass,"n_ice", "number of ice layers",iice)
    1347         CALL put_field(pass,"IR_soil", "Soil IR flux",IRs)
    1348         CALL put_field(pass,"LMO", "Monin-Obukhov Scale",LMO)
    1349         CALL put_field(pass,"surf_water", "Surficial water",rusn)
    1350         CALL put_field(pass,"snow_buffer", "Snow buffer layer",Bufs)
    1351         CALL put_field(pass,"alb_1", "albedo sw",alb1)
    1352         CALL put_field(pass,"alb_2", "albedo nIR",alb2)
    1353         CALL put_field(pass,"alb_3", "albedo fIR",alb3)
    1354         CALL put_field(pass,"to_ice", "Snow passed to ice",toic)
    1355 
    1356 
    1357 
    1358         DO isn = 1,nsno
    1359           IF (isn.LE.99) THEN
    1360             WRITE(str2,'(i2.2)') isn
    1361             CALL put_field(pass,"AGESNOW"//str2, &
    1362                          "Age de la neige layer No."//str2, &
    1363                          agsn(:,isn))
    1364           ELSE
    1365             PRINT*, "Trop de couches"
    1366             CALL abort
    1367           ENDIF
     1268    !======================================================================
     1269    SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat)
     1270
     1271
     1272
     1273        !======================================================================
     1274        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
     1275        ! Objet: Ecriture de l'etat de redemarrage pour SISVAT
     1276        !======================================================================
     1277        USE mod_grid_phy_lmdz
     1278        USE mod_phys_lmdz_para
     1279        USE iostart
     1280        USE VAR_SV
     1281        USE VARxSV
     1282        USE VARySV !hj tmp 12 03 2010
     1283        USE VARtSV
     1284        USE indice_sol_mod
     1285        USE dimphy
     1286
     1287        IMPLICIT none
     1288
     1289        include "netcdf.inc"
     1290        !    include "indicesol.h"
     1291        !    include "dimsoil.h"
     1292        include "clesphys.h"
     1293        include "thermcell.h"
     1294        include "compbl.h"
     1295
     1296        !======================================================================
     1297
     1298        CHARACTER(LEN = *) :: fichnom
     1299        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
     1300        REAL, DIMENSION(klon), INTENT(IN) :: rlon
     1301        REAL, DIMENSION(klon), INTENT(IN) :: rlat
     1302
     1303        ! les variables globales ecrites dans le fichier restart
     1304        REAL, DIMENSION(klon) :: isno
     1305        REAL, DIMENSION(klon) :: ispi
     1306        REAL, DIMENSION(klon) :: iice
     1307        REAL, DIMENSION(klon, nsnowmx) :: isto
     1308
     1309        REAL, DIMENSION(klon, nsismx) :: Tsis
     1310        REAL, DIMENSION(klon, nsismx) :: eta
     1311        REAL, DIMENSION(klon, nsnowmx) :: dzsn
     1312        REAL, DIMENSION(klon, nsismx) :: ro
     1313        REAL, DIMENSION(klon, nsnowmx) :: G1sn
     1314        REAL, DIMENSION(klon, nsnowmx) :: G2sn
     1315        REAL, DIMENSION(klon, nsnowmx) :: agsn
     1316        REAL, DIMENSION(klon) :: IRs
     1317        REAL, DIMENSION(klon) :: LMO
     1318        REAL, DIMENSION(klon) :: rusn
     1319        REAL, DIMENSION(klon) :: toic
     1320        REAL, DIMENSION(klon) :: Bufs
     1321        REAL, DIMENSION(klon) :: alb1, alb2, alb3
     1322
     1323        INTEGER isl, ikl, i, isn, ierr
     1324        CHARACTER (len = 2) :: str2
     1325        INTEGER :: pass
     1326
     1327        isno(:) = 0
     1328        ispi(:) = 0
     1329        iice(:) = 0
     1330        IRs(:) = 0.
     1331        LMO(:) = 0.
     1332        eta(:, :) = 0.
     1333        Tsis(:, :) = 0.
     1334        isto(:, :) = 0
     1335        ro(:, :) = 0.
     1336        G1sn(:, :) = 0.
     1337        G2sn(:, :) = 0.
     1338        dzsn(:, :) = 0.
     1339        agsn(:, :) = 0.
     1340        rusn(:) = 0.
     1341        toic(:) = 0.
     1342        Bufs(:) = 0.
     1343        alb1(:) = 0.
     1344        alb2(:) = 0.
     1345        alb3(:) = 0.
     1346
     1347        !***************************************************************************
     1348        ! Uncompress SISVAT output variables for storage
     1349
     1350        DO  ikl = 1, klon
     1351            i = ikl2i(ikl)
     1352            IF (i > 0) THEN
     1353                isno(i) = 1. * isnoSV(ikl)               ! Nb Snow/Ice Lay.
     1354                ispi(i) = 1. * ispiSV(ikl)               ! Nb Supr.Ice Lay.
     1355                iice(i) = 1. * iiceSV(ikl)               ! Nb      Ice Lay.
     1356
     1357                !        IRs(i)        = IRs_SV(ikl)
     1358                !        LMO(i)        = LMO_SV(ikl)
     1359
     1360                DO isl = -nsol, 0                           !
     1361                    eta(i, nsno + 1 - isl) = eta_SV(ikl, isl)            ! Soil Humidity
     1362                    Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl)            ! Soil Temperature
     1363                    ro(i, nsno + 1 - isl) = ro__SV(ikl, isl)            !        [kg/m3]
     1364                END DO
     1365
     1366                DO  isn = 1, nsno
     1367                    isto(i, isn) = 1. * istoSV(ikl, isn)         ! Snow     History
     1368                    ro(i, isn) = ro__SV(ikl, isn)            !        [kg/m3]
     1369                    eta(i, isn) = eta_SV(ikl, isn)            !        [m3/m3]
     1370                    Tsis(i, isn) = TsisSV(ikl, isn)            !            [K]
     1371                    G1sn(i, isn) = G1snSV(ikl, isn)            ! [-]        [-]
     1372                    G2sn(i, isn) = G2snSV(ikl, isn)            ! [-] [0.0001 m]
     1373                    dzsn(i, isn) = dzsnSV(ikl, isn)            !            [m]
     1374                    agsn(i, isn) = agsnSV(ikl, isn)            !          [day]
     1375                END DO
     1376                rusn(i) = rusnSV(ikl)                  ! Surficial Water
     1377                toic(i) = toicSV(ikl)                  ! to ice
     1378                alb1(i) = alb1sv(ikl)
     1379                alb2(i) = alb2sv(ikl)
     1380                alb3(i) = alb3sv(ikl)
     1381                !        Bufs(i)       = BufsSV(ikl)
     1382            END IF
     1383        END DO
     1384
     1385        CALL open_restartphy(fichnom)
     1386
     1387        DO pass = 1, 2
     1388            CALL put_field(pass, "longitude", &
     1389                    "Longitudes de la grille physique", rlon)
     1390            CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat)
     1391
     1392            CALL put_field(pass, "n_snows", "number of snow/ice layers", isno)
     1393            CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi)
     1394            CALL put_field(pass, "n_ice", "number of ice layers", iice)
     1395            CALL put_field(pass, "IR_soil", "Soil IR flux", IRs)
     1396            CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO)
     1397            CALL put_field(pass, "surf_water", "Surficial water", rusn)
     1398            CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs)
     1399            CALL put_field(pass, "alb_1", "albedo sw", alb1)
     1400            CALL put_field(pass, "alb_2", "albedo nIR", alb2)
     1401            CALL put_field(pass, "alb_3", "albedo fIR", alb3)
     1402            CALL put_field(pass, "to_ice", "Snow passed to ice", toic)
     1403
     1404            DO isn = 1, nsno
     1405                IF (isn.LE.99) THEN
     1406                    WRITE(str2, '(i2.2)') isn
     1407                    CALL put_field(pass, "AGESNOW" // str2, &
     1408                            "Age de la neige layer No." // str2, &
     1409                            agsn(:, isn))
     1410                ELSE
     1411                    PRINT*, "Trop de couches"
     1412                    CALL abort
     1413                ENDIF
     1414            ENDDO
     1415            DO isn = 1, nsno
     1416                IF (isn.LE.99) THEN
     1417                    WRITE(str2, '(i2.2)') isn
     1418                    CALL put_field(pass, "DZSNOW" // str2, &
     1419                            "Snow/ice thickness layer No." // str2, &
     1420                            dzsn(:, isn))
     1421                ELSE
     1422                    PRINT*, "Trop de couches"
     1423                    CALL abort
     1424                ENDIF
     1425            ENDDO
     1426            DO isn = 1, nsno
     1427                IF (isn.LE.99) THEN
     1428                    WRITE(str2, '(i2.2)') isn
     1429                    CALL put_field(pass, "G2SNOW" // str2, &
     1430                            "Snow Property 2, layer No." // str2, &
     1431                            G2sn(:, isn))
     1432                ELSE
     1433                    PRINT*, "Trop de couches"
     1434                    CALL abort
     1435                ENDIF
     1436            ENDDO
     1437            DO isn = 1, nsno
     1438                IF (isn.LE.99) THEN
     1439                    WRITE(str2, '(i2.2)') isn
     1440                    CALL put_field(pass, "G1SNOW" // str2, &
     1441                            "Snow Property 1, layer No." // str2, &
     1442                            G1sn(:, isn))
     1443                ELSE
     1444                    PRINT*, "Trop de couches"
     1445                    CALL abort
     1446                ENDIF
     1447            ENDDO
     1448            DO isn = 1, nsismx
     1449                IF (isn.LE.99) THEN
     1450                    WRITE(str2, '(i2.2)') isn
     1451                    CALL put_field(pass, "ETA" // str2, &
     1452                            "Soil/snow water content layer No." // str2, &
     1453                            eta(:, isn))
     1454                ELSE
     1455                    PRINT*, "Trop de couches"
     1456                    CALL abort
     1457                ENDIF
     1458            ENDDO
     1459            DO isn = 1, nsismx   !nsno
     1460                IF (isn.LE.99) THEN
     1461                    WRITE(str2, '(i2.2)') isn
     1462                    CALL put_field(pass, "RO" // str2, &
     1463                            "Snow density layer No." // str2, &
     1464                            ro(:, isn))
     1465                ELSE
     1466                    PRINT*, "Trop de couches"
     1467                    CALL abort
     1468                ENDIF
     1469            ENDDO
     1470            DO isn = 1, nsismx
     1471                IF (isn.LE.99) THEN
     1472                    WRITE(str2, '(i2.2)') isn
     1473                    CALL put_field(pass, "TSS" // str2, &
     1474                            "Soil/snow temperature layer No." // str2, &
     1475                            Tsis(:, isn))
     1476                ELSE
     1477                    PRINT*, "Trop de couches"
     1478                    CALL abort
     1479                ENDIF
     1480            ENDDO
     1481            DO isn = 1, nsno
     1482                IF (isn.LE.99) THEN
     1483                    WRITE(str2, '(i2.2)') isn
     1484                    CALL put_field(pass, "HISTORY" // str2, &
     1485                            "Snow history layer No." // str2, &
     1486                            isto(:, isn))
     1487                ELSE
     1488                    PRINT*, "Trop de couches"
     1489                    CALL abort
     1490                ENDIF
     1491            ENDDO
     1492
     1493            CALL enddef_restartphy
    13681494        ENDDO
    1369         DO isn = 1,nsno
    1370           IF (isn.LE.99) THEN
    1371             WRITE(str2,'(i2.2)') isn
    1372             CALL put_field(pass,"DZSNOW"//str2, &
    1373                          "Snow/ice thickness layer No."//str2, &
    1374                          dzsn(:,isn))
    1375           ELSE
    1376             PRINT*, "Trop de couches"
    1377             CALL abort
    1378           ENDIF
    1379         ENDDO
    1380         DO isn = 1,nsno
    1381           IF (isn.LE.99) THEN
    1382             WRITE(str2,'(i2.2)') isn
    1383             CALL put_field(pass,"G2SNOW"//str2, &
    1384                          "Snow Property 2, layer No."//str2, &
    1385                          G2sn(:,isn))
    1386           ELSE
    1387             PRINT*, "Trop de couches"
    1388             CALL abort
    1389           ENDIF
    1390         ENDDO
    1391         DO isn = 1,nsno
    1392           IF (isn.LE.99) THEN
    1393             WRITE(str2,'(i2.2)') isn
    1394             CALL put_field(pass,"G1SNOW"//str2, &
    1395                          "Snow Property 1, layer No."//str2, &
    1396                          G1sn(:,isn))
    1397           ELSE
    1398             PRINT*, "Trop de couches"
    1399             CALL abort
    1400           ENDIF
    1401         ENDDO
    1402         DO isn = 1,nsismx
    1403           IF (isn.LE.99) THEN
    1404             WRITE(str2,'(i2.2)') isn
    1405             CALL put_field(pass,"ETA"//str2, &
    1406                          "Soil/snow water content layer No."//str2, &
    1407                          eta(:,isn))
    1408           ELSE
    1409             PRINT*, "Trop de couches"
    1410             CALL abort
    1411           ENDIF
    1412         ENDDO
    1413         DO isn = 1,nsismx   !nsno
    1414           IF (isn.LE.99) THEN
    1415             WRITE(str2,'(i2.2)') isn
    1416             CALL put_field(pass,"RO"//str2, &
    1417                            "Snow density layer No."//str2, &
    1418                            ro(:,isn))
    1419           ELSE
    1420             PRINT*, "Trop de couches"
    1421             CALL abort
    1422           ENDIF
    1423         ENDDO
    1424         DO isn = 1,nsismx
    1425           IF (isn.LE.99) THEN
    1426             WRITE(str2,'(i2.2)') isn
    1427             CALL put_field(pass,"TSS"//str2, &
    1428                            "Soil/snow temperature layer No."//str2, &
    1429                            Tsis(:,isn))
    1430           ELSE
    1431             PRINT*, "Trop de couches"
    1432             CALL abort
    1433           ENDIF
    1434         ENDDO
    1435         DO isn = 1,nsno
    1436           IF (isn.LE.99) THEN
    1437             WRITE(str2,'(i2.2)') isn
    1438             CALL put_field(pass,"HISTORY"//str2, &
    1439                            "Snow history layer No."//str2, &
    1440                            isto(:,isn))
    1441           ELSE
    1442             PRINT*, "Trop de couches"
    1443             CALL abort
    1444           ENDIF
    1445         ENDDO
    1446 
    1447       CALL enddef_restartphy
    1448       ENDDO
    1449       CALL close_restartphy
    1450 
    1451 
    1452   END SUBROUTINE sisvatredem
     1495        CALL close_restartphy
     1496
     1497    END SUBROUTINE sisvatredem
    14531498
    14541499END MODULE surf_inlandsis_mod
Note: See TracChangeset for help on using the changeset viewer.