New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14572 – NEMO

Changeset 14572


Ignore:
Timestamp:
2021-03-03T14:08:29+01:00 (3 years ago)
Author:
antsia
Message:

Merge with dev_isf_flx_UKESM_r9321@10026

Location:
branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/CONFIG/SHARED/field_def.xml

    r9163 r14572  
    240240         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature)"   unit="m/s"      /> 
    241241         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)"      unit="m/s"      /> 
    242          <field id="stbl"         long_name="salinity in the Losh tbl"                      unit="0.001"     /> 
    243          <field id="ttbl"         long_name="temperature in the Losh tbl"                   unit="degree_C"     /> 
     242         <field id="stbl"         long_name="salinity in the Losh tbl"                      unit="0.001"    /> 
     243         <field id="ttbl"         long_name="temperature in the Losh tbl"                   unit="degree_C" /> 
     244         <field id="utbl"         long_name="zonal current in the Losh tbl at T point"      unit="m/s"      /> 
     245         <field id="vtbl"         long_name="meridional current in the Losh tbl at T point" unit="m/s"      /> 
     246         <field id="isftfrz"      long_name="freezing point temperature"                    unit="degree_C" /> 
     247         <field id="isfsfrz"      long_name="salinity at the ice-shelf/ocean interface"     unit="0.001"    /> 
     248 
    244249 
    245250         <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 
    … …  
    442447         <field id="ahu_bbl"      long_name="BBL diffusive flux along i-axis"   unit="m3/s" /> 
    443448 
    444          <!-- variable for ice shelves --> 
    445          <field id="utbl"         long_name="zonal current in the Losh tbl"     unit="m/s" /> 
    446  
    447449         <!-- variables available with key_diaar5 --> 
    448450         <field id="u_masstr"     long_name="Ocean Mass X Transport"    standard_name="ocean_mass_x_transport"                          unit="kg/s"        grid_ref="grid_U_3D" /> 
    … …  
    484486         <field id="voce_bbl"     long_name="BBL ocean current along j-axis"    unit="m/s"  /> 
    485487         <field id="ahv_bbl"      long_name="BBL diffusive flux along j-axis"   unit="m3/s" /> 
    486  
    487          <!-- variable for ice shelves --> 
    488          <field id="vtbl"         long_name="meridional current in the Losh tbl"   unit="m/s" /> 
    489488 
    490489         <!-- variables available with key_diaar5 --> 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/CONFIG/SHARED/namelist_ref

    r8447 r14572  
    471471!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! 
    472472! nn_isf == 4 
    473    sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
    474473   sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
    475474! nn_isf == 3 
    … …  
    480479! nn_isf == 2 
    481480   sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,   '' 
    482 ! for all case 
    483    ln_divisf   = .true.  ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?) 
    484481! only for nn_isf = 1 or 2 
    485482   rn_gammat0  = 1.0e-4   ! gammat coefficient used in blk formula 
    … …  
    489486   rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008) 
    490487                          ! 0 => thickness of the tbl = thickness of the first wet cell 
    491    ln_conserve = .true.   ! conservative case (take into account meltwater advection) 
    492488   nn_gammablk = 1        ! 0 = cst Gammat (= gammat/s) 
    493489                          ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r8427 r14572  
    259259   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: first wet T-, U-, V-, F- ocean level (ISF) 
    260260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                       (ISF) 
    261    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask                   !: surface domain T-point mask  
     261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssisfmask       !: surface domain T-point mask  
    262262 
    263263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    … …  
    403403      ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),     & 
    404404         &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
    405          &     mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 
     405         &     mikf(jpi,jpj), ssmask(jpi,jpj), ssisfmask(jpi,jpj), STAT=ierr(10) ) 
    406406 
    407407      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r8280 r14572  
    218218      ENDIF 
    219219#endif 
     220 
     221      ! in case of iceshelf definition of iceshelf surface mask (1 under 
     222      ! iceshelf and 0 elsewhere 
     223      ssisfmask = (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     224 
    220225!!gm end 
    221226 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r9321 r14572  
    197197       
    198198      ! note that mbkt is set to 1 over land ==> use surface tmask 
    199       zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp ) 
     199      zprt(:,:) = ssisfmask(:,:) 
     200      CALL iom_rstput( 0, 0, inum4, 'maskisf', zprt, ktype = jp_i1 )       !    ! surface ice shelf mask  
     201      zprt(:,:) = ssmask(:,:) * mbkt(:,:) 
    200202      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points 
    201203      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp ) 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r6487 r14572  
    229229 
    230230      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field) 
    231       IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
     231      IF( nn_isf /= 0 )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
    232232      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence) 
    233233       
    … …  
    328328 
    329329      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field) 
    330       IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
     330      IF( nn_isf /=  0)   CALL sbc_isf_div( hdivn )                            ! ice shelf (update hdivn field) 
    331331      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    332332      ! 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r8046 r14572  
    4747   LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice) 
    4848   INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3) 
    49    INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4)  
     49   INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4) 
    5050   INTEGER , PUBLIC ::   nn_ice_embd    !: flag for levitating/embedding sea-ice in the ocean 
    5151   !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect) 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r8046 r14572  
    1818   USE eosbn2          ! equation of state 
    1919   USE sbc_oce         ! surface boundary condition: ocean fields 
     20   USE zdfbfr          ! 
     21   ! 
     22   USE in_out_manager  ! I/O manager 
     23   USE iom             ! I/O manager library 
     24   USE fldread         ! read input field at current time step 
    2025   USE lbclnk          ! 
    21    USE iom             ! I/O manager library 
    22    USE in_out_manager  ! I/O manager 
    2326   USE wrk_nemo        ! Memory allocation 
    2427   USE timing          ! Timing 
    2528   USE lib_fortran     ! glob_sum 
    26    USE zdfbfr 
    27    USE fldread         ! read input field at current time step 
    28    USE lib_fortran, ONLY: glob_sum 
    2929 
    3030   IMPLICIT NONE 
    … …  
    3535   ! public in order to be able to output then  
    3636 
    37    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
    38    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf 
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc  !: before and now T & S isf contents [K.m/s & PSU.m/s]  
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf, qhcisf          !: latent heat and heat content flux from ice shelf  [W/m2] 
    3939   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m] 
    40    LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence  
    41    INTEGER , PUBLIC ::   nn_isfblk                   !:  
    42    INTEGER , PUBLIC ::   nn_gammablk                 !: 
    43    LOGICAL , PUBLIC ::   ln_conserve                 !: 
    44    REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient 
    45    REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient  
    46    REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence 
     40   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting 
     41   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed 
     42   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient [] 
     43   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient [] 
    4744 
    4845   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3 
    49    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl 
     46   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl  [m] 
    5047   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/thickness of tbl 
    5148   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by tbl  
    5249   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2 
    5350   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 
    54    INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base 
    55  
    56  
    57    REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ? 
    58    REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ? 
    59    REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ? 
    60    REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ? 
    61    REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ? 
     51   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)      ::  misfkt, misfkb         !:Level of ice shelf base 
     52 
     53   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! specific heat of ice shelf             [J/kg/K] 
     54   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! heat diffusivity through the ice-shelf [m2/s] 
     55   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! volumic mass of ice shelf              [kg/m3] 
     56   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! air temperature on top of ice shelf    [C] 
     57   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! latent heat of fusion of ice shelf     [J/kg] 
    6258 
    6359!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 
    64    CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files 
    65    TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read 
    66    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf 
    67    TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read 
    68    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf            
    69    TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read 
     60   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files 
     61   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read 
     62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 
     63   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read 
     64   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf            
     65   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read 
     66   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read 
     67   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read 
    7068    
    7169   !! * Substitutions 
    … …  
    7674   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    7775   !!---------------------------------------------------------------------- 
    78  
    7976CONTAINS 
    8077  
    8178  SUBROUTINE sbc_isf(kt) 
    82     INTEGER, INTENT(in)          ::   kt         ! ocean time step 
    83     INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror 
    84     INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer 
    85     REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum 
    86     REAL(wp)                     ::   rmin 
    87     REAL(wp)                     ::   zhk 
    88     REAL(wp)                     ::   zt_frz, zpress 
    89     CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file 
    90     CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file 
    91     CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale 
    92     INTEGER           ::   ios           ! Local integer output status for namelist read 
    93  
    94     REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 
    95     REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d 
    96       ! 
    9779      !!--------------------------------------------------------------------- 
    98       NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, & 
    99                          & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
     80      !!                  ***  ROUTINE sbc_isf  *** 
     81      !! 
     82      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf  
     83      !!              melting and freezing  
     84      !! 
     85      !! ** Method  :  4 parameterizations are available according to nn_isf  
     86      !!               nn_isf = 1 : Realistic ice_shelf formulation 
     87      !!                        2 : Beckmann & Goose parameterization 
     88      !!                        3 : Specified runoff in deptht (Mathiot & al. ) 
     89      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
     90      !!---------------------------------------------------------------------- 
     91      INTEGER, INTENT(in)          ::   kt         ! ocean time step 
     92      INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror 
     93      INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer 
     94      REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum 
     95      REAL(wp)                     ::   rmin 
     96      REAL(wp)                     ::   zhk 
     97      CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file 
     98      CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file 
     99      CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale 
     100      INTEGER                      :: ios                         ! Local integer output status for namelist read 
     101 
     102      REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 
     103      REAL(wp), DIMENSION(:,:  ), POINTER :: ztfrz, zdep         ! freezing temperature (ztfrz) at depth (zdep)  
     104      ! 
     105      !!--------------------------------------------------------------------- 
     106      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, & 
     107                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
    100108      ! 
    101109      ! 
    … …  
    121129         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk 
    122130         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl 
    123          IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf  
    124131         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
    125132         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2  
    126          IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM) 
    127             rdivisf = 1._wp 
    128          ELSE 
    129             rdivisf = 0._wp 
    130          END IF 
    131133         ! 
    132134         ! Allocate public variable 
    … …  
    185187             
    186188            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
    187             ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 
     189            ALLOCATE( sf_fwfisf(1), STAT=ierror ) 
    188190            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
    189             ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 
    190191            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    191             !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' ) 
    192192         END IF 
    193193          
    … …  
    206206 
    207207      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 
     208         ! allocation 
     209         CALL wrk_alloc( jpi,jpj, ztfrz, zdep  ) 
    208210 
    209211         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 
    … …  
    229231 
    230232         ! compute salf and heat flux 
    231          IF (nn_isf == 1) THEN 
    232             ! realistic ice shelf formulation 
     233         SELECT CASE ( nn_isf ) 
     234         CASE ( 1 )    ! realistic ice shelf formulation 
    233235            ! compute T/S/U/V for the top boundary layer 
    234236            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
    235237            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 
    236             CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U') 
    237             CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V') 
     238            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U') 
     239            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V') 
    238240            ! iom print 
    239241            CALL iom_put('ttbl',ttbl(:,:)) 
    … …  
    244246            CALL sbc_isf_cav (kt) 
    245247 
    246          ELSE IF (nn_isf == 2) THEN 
    247             ! Beckmann and Goosse parametrisation  
     248         CASE ( 2 )    ! Beckmann and Goosse parametrisation  
    248249            stbl(:,:)   = soce 
    249250            CALL sbc_isf_bg03(kt) 
    250251 
    251          ELSE IF (nn_isf == 3) THEN 
    252             ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
     252         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    253253            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    254             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
     254            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf flux from the isf (fwfisf <0 mean melting)  
    255255 
    256256            IF( lk_oasis) THEN 
    … …  
    294294            ENDIF 
    295295 
    296             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
     296            qisf(:,:)   = fwfisf(:,:) * rlfusisf              ! heat        flux 
    297297            stbl(:,:)   = soce 
    298298 
    299          ELSE IF (nn_isf == 4) THEN 
    300             ! specified fwf and heat flux forcing beneath the ice shelf 
     299         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf 
    301300            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   ) 
    302             !CALL fld_read ( kt, nn_fsbc, sf_qisf   ) 
    303             fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
     301            fwfisf(:,:) =   sf_fwfisf(1)%fnow(:,:,1)            ! fwf  flux from the isf (fwfisf <0 mean melting) 
    304302 
    305303            IF( lk_oasis) THEN 
    … …  
    343341            ENDIF 
    344342 
    345             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
    346             !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux 
     343            qisf(:,:)   = fwfisf(:,:) * rlfusisf              ! heat        flux 
    347344            stbl(:,:)   = soce 
    348345 
    349          END IF 
     346         END SELECT 
     347 
     348         ! compute heat content 
     349         ! heat content flux for nn_isf==1 is done in sbcisf_cav 
     350         SELECT CASE ( nn_isf ) 
     351         CASE ( 2, 3, 4 ) 
     352            DO jj = 1,jpj 
     353               DO ji = 1,jpi 
     354                  zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj)) 
     355               END DO 
     356            END DO 
     357            CALL eos_fzp( stbl(:,:), ztfrz(:,:), zdep(:,:) ) 
     358            qhcisf(:,:) = - fwfisf(:,:) * ztfrz(:,:) * rcp 
     359         END SELECT 
     360          
    350361         ! compute tsc due to isf 
    351          ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable). 
    352 !         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
    353          zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 
    354          risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 ! 
    355           
    356          ! salt effect already take into account in vertical advection 
    357          risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 
     362         risf_tsc(:,:,jp_tem) = (qisf(:,:) + qhcisf(:,:)) * r1_rau0_rcp 
     363         risf_tsc(:,:,jp_sal) = 0.0_wp 
    358364 
    359365         ! output 
    360          IF( iom_use('qlatisf' ) )   CALL iom_put('qlatisf', qisf) 
    361          IF( iom_use('fwfisf'  ) )   CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce ) 
    362  
    363          ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now 
    364          fwfisf(:,:) = rdivisf * fwfisf(:,:)          
    365   
     366         IF( iom_use('qlatisf' ) )   CALL iom_put('qlatisf', qisf  ) 
     367         IF( iom_use('fwfisf'  ) )   CALL iom_put('fwfisf' , fwfisf) 
     368         IF( iom_use('qhcisf'  ) )   CALL iom_put('qhcisf' , qhcisf) 
     369 
    366370         ! lbclnk 
    367371         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 
    368372         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 
    369          CALL lbc_lnk(fwfisf(:,:)   ,'T',1.) 
    370          CALL lbc_lnk(qisf(:,:)     ,'T',1.) 
     373         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.) 
     374         CALL lbc_lnk(qisf(:,:)           ,'T',1.) 
    371375 
    372376!============================================================================================================================================= 
    373          IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 
     377         IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d')) THEN 
    374378            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 
    375             CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        ) 
    376379 
    377380            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s) 
    378381            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2) 
    379382            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2) 
    380             zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2) 
    381383 
    382384            DO jj = 1,jpj 
    … …  
    385387                  ikb = misfkb(ji,jj) 
    386388                  DO jk = ikt, ikb - 1 
    387                      zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
    388                      zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
    389                      zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
     389                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
     390                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + qhcisf(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
     391                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf  (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
    390392                  END DO 
    391                   zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
    392                   zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
    393                   zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
     393                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
     394                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + qhcisf(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
     395                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf  (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
    394396               END DO 
    395397            END DO 
    … …  
    398400            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 
    399401            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 
    400             CALL iom_put('qhcisf'   , zqhcisf2d (:,:  )) 
    401402 
    402403            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 
    403             CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        ) 
    404404         END IF 
    405405!============================================================================================================================================= 
    406406 
    407          IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
     407         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    ! 
    408408            IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
    409409                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 
    … …  
    416416               risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 
    417417            END IF 
    418          ENDIF 
     418         END IF 
    419419         !  
     420         ! deallocation 
     421         CALL wrk_dealloc( jpi,jpj, ztfrz, zdep  ) 
    420422      END IF 
    421    
     423      !   
    422424  END SUBROUTINE sbc_isf 
     425 
    423426 
    424427  INTEGER FUNCTION sbc_isf_alloc() 
    … …  
    433436               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), & 
    434437               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , & 
    435                &    STAT= sbc_isf_alloc ) 
     438               &    qhcisf(jpi,jpj), STAT= sbc_isf_alloc ) 
    436439         ! 
    437          IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc ) 
     440         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc ) 
    438441         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 
    439442         ! 
    440       ENDIF 
     443      END IF 
    441444  END FUNCTION 
    442445 
    443446  SUBROUTINE sbc_isf_bg03(kt) 
    444    !!========================================================================== 
    445    !!                 *** SUBROUTINE sbcisf_bg03  *** 
    446    !! add net heat and fresh water flux from ice shelf melting 
    447    !! into the adjacent ocean using the parameterisation by 
    448    !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
    449    !!     interaction for climate models", Ocean Modelling 5(2003) 157-170. 
    450    !!  (hereafter BG) 
    451    !!========================================================================== 
    452    !!---------------------------------------------------------------------- 
    453    !!   sbc_isf_bg03      : routine called from sbcmod 
    454    !!---------------------------------------------------------------------- 
    455    !! 
    456    !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting 
    457    !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling 
    458    !! 
    459    !! History : 
    460    !!      !  06-02  (C. Wang) Original code 
    461    !!---------------------------------------------------------------------- 
    462  
    463     INTEGER, INTENT ( in ) :: kt 
    464  
    465     INTEGER :: ji, jj, jk, jish  !temporary integer 
    466     INTEGER :: ijkmin 
    467     INTEGER :: ii, ij, ik  
    468     INTEGER :: inum 
    469  
    470     REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m 
    471     REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m 
    472     REAL(wp) :: zt_frz      ! freezing point temperature at depth z 
    473     REAL(wp) :: zpress      ! pressure to compute the freezing point in depth 
    474      
    475     !!---------------------------------------------------------------------- 
    476     IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
    477      ! 
    478  
    479     ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run ) 
    480     DO ji = 1, jpi 
    481        DO jj = 1, jpj 
    482           ik = misfkt(ji,jj) 
    483           !! Initialize arrays to 0 (each step) 
    484           zt_sum = 0.e0_wp 
    485           IF ( ik .GT. 1 ) THEN 
    486     ! 3. -----------the average temperature between 200m and 600m --------------------- 
    487              DO jk = misfkt(ji,jj),misfkb(ji,jj) 
    488              ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 
    489              ! after verif with UNESCO, wrong sign in BG eq. 2 
    490              ! Calculate freezing temperature 
    491                 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04  
    492                 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress)  
    493                 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp 
    494              ENDDO 
    495              zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
    496      
    497     ! 4. ------------Net heat flux and fresh water flux due to the ice shelf 
    498           ! For those corresponding to zonal boundary     
    499              qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
    500                          & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik)  
     447      !!--------------------------------------------------------------------- 
     448      !!                  ***  ROUTINE sbc_isf_bg03  *** 
     449      !! 
     450      !! ** Purpose : add net heat and fresh water flux from ice shelf melting 
     451      !!          into the adjacent ocean 
     452      !! 
     453      !! ** Method  :   See reference 
     454      !! 
     455      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
     456      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170. 
     457      !!         (hereafter BG) 
     458      !! History : 
     459      !!         06-02  (C. Wang) Original code 
     460      !!---------------------------------------------------------------------- 
     461      INTEGER, INTENT ( in ) :: kt 
     462      ! 
     463      INTEGER  :: ji, jj, jk ! dummy loop index 
     464      INTEGER  :: ik         ! current level 
     465      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m 
     466      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m 
     467      REAL(wp) :: ztfrz     ! freezing point temperature at depth z 
     468      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth 
     469      !!---------------------------------------------------------------------- 
     470 
     471      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
     472      ! 
     473      DO ji = 1, jpi 
     474         DO jj = 1, jpj 
     475            ik = misfkt(ji,jj) 
     476            !! Initialize arrays to 0 (each step) 
     477            zt_sum = 0.e0_wp 
     478            IF ( ik > 1 ) THEN 
     479               ! 1. -----------the average temperature between 200m and 600m --------------------- 
     480               DO jk = misfkt(ji,jj),misfkb(ji,jj) 
     481                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 
     482                  ! after verif with UNESCO, wrong sign in BG eq. 2 
     483                  ! Calculate freezing temperature 
     484                  CALL eos_fzp(stbl(ji,jj), ztfrz, gdept_n(ji,jj,jk))  
     485                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-ztfrz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
     486               END DO 
     487               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
     488               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 
     489               ! For those corresponding to zonal boundary     
     490               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
     491                           & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk) 
    501492              
    502              fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                   
    503              fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
    504              !add to salinity trend 
    505           ELSE 
    506              qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 
    507           END IF 
    508        ENDDO 
    509     ENDDO 
    510     ! 
    511     IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03') 
     493               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                   
     494               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
     495               !add to salinity trend 
     496            ELSE 
     497               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 
     498            END IF 
     499         END DO 
     500      END DO 
     501      ! 
     502      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03') 
     503      ! 
    512504  END SUBROUTINE sbc_isf_bg03 
    513505 
    … …  
    527519      INTEGER, INTENT(in)          ::   kt         ! ocean time step 
    528520      ! 
    529       LOGICAL :: ln_isomip = .true. 
    530       REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti 
    531       REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d  
    532       !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf 
     521      INTEGER  ::   ji, jj     ! dummy loop indices 
     522      INTEGER  ::   nit 
    533523      REAL(wp) ::   zlamb1, zlamb2, zlamb3 
    534524      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 
    535       REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 
    536       REAL(wp) ::   zfwflx, zhtflx, zhtflx_b 
    537       REAL(wp) ::   zgammat, zgammas 
    538       REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==! 
    539       INTEGER  ::   ji, jj     ! dummy loop indices 
    540       INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers 
    541       INTEGER  ::   ierror     ! return error code 
    542       LOGICAL  ::   lit=.TRUE. 
    543       INTEGER  ::   nit 
     525      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zcfac 
     526      REAL(wp) ::   zeps = 1.e-20_wp         
     527      REAL(wp) ::   zerr 
     528      REAL(wp), DIMENSION(:,:), POINTER ::   ztfrz, zsfrz 
     529      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas  
     530      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b 
     531      LOGICAL  ::   lit 
    544532      !!--------------------------------------------------------------------- 
    545       ! 
    546       ! coeficient for linearisation of tfreez 
    547       zlamb1=-0.0575 
    548       zlamb2=0.0901 
    549       zlamb3=-7.61e-04 
     533      ! coeficient for linearisation of potential tfreez 
     534      ! Crude approximation for pressure (but commonly used) 
     535      IF ( ln_useCT ) THEN   ! linearisation from Jourdain et al. (2017) 
     536         zlamb1 =-0.0564_wp 
     537         zlamb2 = 0.0773_wp 
     538         zlamb3 =-7.8633e-8 * grav * rau0 
     539      ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015) 
     540         zlamb1 =-0.0573_wp 
     541         zlamb2 = 0.0832_wp 
     542         zlamb3 =-7.53e-8 * grav * rau0 
     543      ENDIF 
     544      ! 
    550545      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav') 
    551546      ! 
    552       CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 
    553  
    554       zcfac=0.0_wp  
    555       IF (ln_conserve)  zcfac=1.0_wp 
    556       zpress(:,:)=0.0_wp 
    557       zgammat2d(:,:)=0.0_wp 
    558       zgammas2d(:,:)=0.0_wp 
    559       ! 
    560       ! 
    561 !CDIR COLLAPSE 
    562       DO jj = 1, jpj 
    563          DO ji = 1, jpi 
    564             ! Crude approximation for pressure (but commonly used) 
    565             ! 1e-04 to convert from Pa to dBar 
    566             zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04 
     547      CALL wrk_alloc( jpi,jpj, ztfrz  , zsfrz  , zgammat, zgammas  ) 
     548      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 
     549 
     550      ! initialisation 
     551      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 
     552      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp     
     553      zfwflx (:,:) = 0.0_wp     ; zsfrz(:,:) = 0.0_wp 
     554 
     555      ! compute ice shelf melting 
     556      nit = 1 ; lit = .TRUE. 
     557      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine 
     558         SELECT CASE ( nn_isfblk ) 
     559         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 
     560            ! Calculate freezing temperature 
     561            CALL eos_fzp( stbl(:,:), ztfrz(:,:), risfdep(:,:) ) 
     562 
     563            ! compute gammat every where (2d) 
     564            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     565             
     566            ! compute upward heat flux zhtflx and upward water flux zwflx 
     567            DO jj = 1, jpj 
     568               DO ji = 1, jpi 
     569                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-ztfrz(ji,jj)) 
     570                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 
     571               END DO 
     572            END DO 
     573 
     574            ! Compute heat flux and upward fresh water flux 
     575            qisf  (:,:) = - zhtflx(:,:) * ssisfmask(:,:) 
     576            fwfisf(:,:) =   zfwflx(:,:) * ssisfmask(:,:) 
     577 
     578         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
     579            ! compute gammat every where (2d) 
     580            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     581 
     582            ! compute upward heat flux zhtflx and upward water flux zwflx 
     583            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 
     584            DO jj = 1, jpj 
     585               DO ji = 1, jpi 
     586                  ! compute coeficient to solve the 2nd order equation 
     587                  zeps1 = rcp*rau0*zgammat(ji,jj) 
     588                  zeps2 = rlfusisf*rau0*zgammas(ji,jj) 
     589                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 
     590                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 
     591                  zeps6 = zeps4-ttbl(ji,jj) 
     592                  zeps7 = zeps4-tsurf 
     593                  zaqe  = zlamb1 * (zeps1 + zeps3) 
     594                  zaqer = 0.5_wp/MIN(zaqe,-zeps) 
     595                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2 
     596                  zcqe  = zeps2*stbl(ji,jj) 
     597                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe                
     598 
     599                  ! Presumably zdis can never be negative because gammas is very small compared to gammat 
     600                  ! compute s freeze 
     601                  zsfrz(ji,jj)=(-zbqe-SQRT(zdis))*zaqer 
     602                  IF ( zsfrz(ji,jj) < 0.0_wp ) zsfrz(ji,jj)=(-zbqe+SQRT(zdis))*zaqer 
     603 
     604                  ! compute t freeze (eq. 22) 
     605                  ztfrz(ji,jj)=( zeps4+zlamb1*zsfrz(ji,jj) ) * ssisfmask(ji,jj) 
     606   
     607                  ! zfwflx is upward water flux 
     608                  ! zhtflx is upward heat flux (out of ocean) 
     609                  ! compute the upward water and heat flux (eq. 28 and eq. 29) 
     610                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz(ji,jj)-stbl(ji,jj)) / MAX(zsfrz(ji,jj),zeps) 
     611                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - ztfrz(ji,jj) )  
     612               END DO 
     613            END DO 
    567614            ! 
    568          END DO 
     615            ! mask heat and water flux 
     616            fwfisf(:,:) =   zfwflx(:,:) * ssisfmask(:,:) 
     617            qisf  (:,:) = - zhtflx(:,:) * ssisfmask(:,:) 
     618            ! 
     619            ! compute heat content flux 
     620            qhcisf(:,:) = - fwfisf(:,:) * ztfrz(:,:) * rcp  
     621 
     622         END SELECT 
     623 
     624         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 
     625         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE. 
     626         ELSE                            
     627            ! check total number of iteration 
     628            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
     629            ELSE                 ; nit = nit + 1 
     630            END IF 
     631 
     632            ! compute error between 2 iterations 
     633            ! if needed save gammat and compute zhtflx_b for next iteration 
     634            zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 
     635            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 
     636            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:) 
     637            END IF 
     638         END IF 
    569639      END DO 
    570  
    571 ! Calculate in-situ temperature (ref to surface) 
    572       zti(:,:)=tinsitu( ttbl, stbl, zpress ) 
    573 ! Calculate freezing temperature 
    574       CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress ) 
    575  
    576        
    577       zhtflx=0._wp ; zfwflx=0._wp 
    578       IF (nn_isfblk == 1) THEN 
    579          DO jj = 1, jpj 
    580             DO ji = 1, jpi 
    581                IF (mikt(ji,jj) > 1 ) THEN 
    582                   nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp 
    583                   DO WHILE ( lit ) 
    584 ! compute gamma 
    585                      CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 
    586 ! zhtflx is upward heat flux (out of ocean) 
    587                      zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj)) 
    588 ! zwflx is upward water flux 
    589                      zfwflx = - zhtflx/lfusisf 
    590 ! test convergence and compute gammat 
    591                      IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 
    592  
    593                      nit = nit + 1 
    594                      IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    595  
    596 ! save gammat and compute zhtflx_b 
    597                      zgammat2d(ji,jj)=zgammat 
    598                      zhtflx_b = zhtflx 
    599                   END DO 
    600  
    601                   qisf(ji,jj) = - zhtflx 
    602 ! For genuine ISOMIP protocol this should probably be something like 
    603                   fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps)) 
    604                ELSE 
    605                   fwfisf(ji,jj) = 0._wp 
    606                   qisf(ji,jj)   = 0._wp 
    607                END IF 
    608             ! 
    609             END DO 
    610          END DO 
    611  
    612       ELSE IF (nn_isfblk == 2 ) THEN 
    613  
    614 ! More complicated 3 equation thermodynamics as in MITgcm 
    615 !CDIR COLLAPSE 
    616          DO jj = 2, jpj 
    617             DO ji = 2, jpi 
    618                IF (mikt(ji,jj) > 1 ) THEN 
    619                   nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp 
    620                   DO WHILE ( lit ) 
    621                      CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 
    622  
    623                      zeps1=rcp*rau0*zgammat 
    624                      zeps2=lfusisf*rau0*zgammas 
    625                      zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj) 
    626                      zeps4=zlamb2+zlamb3*risfdep(ji,jj) 
    627                      zeps6=zeps4-zti(ji,jj) 
    628                      zeps7=zeps4-tsurf 
    629                      zaqe=zlamb1 * (zeps1 + zeps3) 
    630                      zaqer=0.5/zaqe 
    631                      zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 
    632                      zcqe=zeps2*stbl(ji,jj) 
    633                      zdis=zbqe*zbqe-4.0*zaqe*zcqe                
    634 ! Presumably zdis can never be negative because gammas is very small compared to gammat 
    635                      zsfrz=(-zbqe-SQRT(zdis))*zaqer 
    636                      IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
    637                      zfrz(ji,jj)=zeps4+zlamb1*zsfrz 
    638    
    639 ! zfwflx is upward water flux 
    640                      zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz ) 
    641 ! zhtflx is upward heat flux (out of ocean) 
    642 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 
    643                      zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) )  
    644 ! zwflx is upward water flux 
    645 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 
    646                      zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 
    647 ! test convergence and compute gammat 
    648                      IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 
    649  
    650                      nit = nit + 1 
    651                      IF (nit .GE. 51) THEN 
    652                         WRITE(numout,*) "sbcisf : too many iteration ... ", & 
    653                             &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea 
    654                         CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    655                      END IF 
    656 ! save gammat and compute zhtflx_b 
    657                      zgammat2d(ji,jj)=zgammat 
    658                      zgammas2d(ji,jj)=zgammas 
    659                      zhtflx_b = zhtflx 
    660  
    661                   END DO 
    662 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 
    663                   qisf(ji,jj) = - zhtflx  
    664 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 
    665                   fwfisf(ji,jj) = zfwflx  
    666                ELSE 
    667                   fwfisf(ji,jj) = 0._wp 
    668                   qisf(ji,jj)   = 0._wp 
    669                ENDIF 
    670                ! 
    671             END DO 
    672          END DO 
    673       ENDIF 
    674       ! lbclnk 
    675       CALL lbc_lnk(zgammas2d(:,:),'T',1.) 
    676       CALL lbc_lnk(zgammat2d(:,:),'T',1.) 
    677       ! output 
    678       CALL iom_put('isfgammat', zgammat2d) 
    679       CALL iom_put('isfgammas', zgammas2d) 
    680       ! 
    681       CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 
     640      ! 
     641      CALL iom_put('isftfrz'  , ztfrz  ) 
     642      CALL iom_put('isfsfrz'  , zsfrz  ) 
     643      CALL iom_put('isfgammat', zgammat) 
     644      CALL iom_put('isfgammas', zgammas) 
     645      !  
     646      CALL wrk_dealloc( jpi,jpj, ztfrz  , zgammat, zgammas  ) 
     647      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 
    682648      ! 
    683649      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav') 
    684  
     650      ! 
    685651   END SUBROUTINE sbc_isf_cav 
    686652 
    687    SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit ) 
     653   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 
    688654      !!---------------------------------------------------------------------- 
    689655      !! ** Purpose    : compute the coefficient echange for heat flux 
    … …  
    694660      !!                Jenkins et al., 2010, JPO, p2298-2312 
    695661      !!--------------------------------------------------------------------- 
    696       REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf 
    697       INTEGER , INTENT(in)    :: ji,jj 
    698       LOGICAL , INTENT(inout) :: lit 
    699  
    700       INTEGER  :: ikt                 ! loop index 
    701       REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity 
     662      REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 
     663      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 
     664      ! 
     665      INTEGER  :: ikt                         
     666      INTEGER  :: ji, jj                     ! loop index 
     667      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity 
    702668      REAL(wp) :: zdku, zdkv                 ! U, V shear  
    703669      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number  
    … …  
    709675      REAL(wp) :: zcoef                      ! temporary coef 
    710676      REAL(wp) :: zdep 
    711       REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant 
    712       REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number 
    713       REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1) 
    714       REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s 
     677      REAL(wp) :: zeps = 1.0e-20_wp     
     678      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant 
     679      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 
    715680      REAL(wp), DIMENSION(2) :: zts, zab 
    716681      !!--------------------------------------------------------------------- 
    717       ! 
    718       IF( nn_gammablk == 0 ) THEN 
    719       !! gamma is constant (specified in namelist) 
    720          gt = rn_gammat0 
    721          gs = rn_gammas0 
    722          lit = .FALSE. 
    723       ELSE IF ( nn_gammablk == 1 ) THEN 
    724       !! gamma is assume to be proportional to u*  
    725       !! WARNING in case of Losh 2008 tbl parametrization,  
    726       !! you have to used the mean value of u in the boundary layer)  
    727       !! not yet coded 
    728       !! Jenkins et al., 2010, JPO, p2298-2312 
    729          ikt = mikt(ji,jj) 
    730       !! Compute U and V at T points 
    731    !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) ) 
    732    !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) ) 
    733           zut = utbl(ji,jj) 
    734           zvt = vtbl(ji,jj) 
    735  
    736       !! compute ustar 
    737          zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 
    738       !! Compute mean value over the TBL 
    739  
    740       !! Compute gammats 
    741          gt = zustar * rn_gammat0 
    742          gs = zustar * rn_gammas0 
    743          lit = .FALSE. 
    744       ELSE IF ( nn_gammablk == 2 ) THEN 
    745       !! gamma depends of stability of boundary layer 
    746       !! WARNING in case of Losh 2008 tbl parametrization,  
    747       !! you have to used the mean value of u in the boundary layer)  
    748       !! not yet coded 
    749       !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
    750       !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
     682      CALL wrk_alloc( jpi,jpj, zustar ) 
     683      ! 
     684      SELECT CASE ( nn_gammablk ) 
     685      CASE ( 0 ) ! gamma is constant (specified in namelist) 
     686         !! ISOMIP formulation (Hunter et al, 2006) 
     687         pgt(:,:) = rn_gammat0 
     688         pgs(:,:) = rn_gammas0 
     689 
     690      CASE ( 1 ) ! gamma is assume to be proportional to u* 
     691         !! Jenkins et al., 2010, JPO, p2298-2312 
     692         !! Adopted by Asay-Davis et al. (2015) 
     693 
     694         !! compute ustar (eq. 24) 
     695         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
     696 
     697         !! Compute gammats 
     698         pgt(:,:) = zustar(:,:) * rn_gammat0 
     699         pgs(:,:) = zustar(:,:) * rn_gammas0 
     700       
     701      CASE ( 2 ) ! gamma depends of stability of boundary layer 
     702         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
     703         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
     704         !! compute ustar 
     705         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
     706 
     707         !! compute Pr and Sc number (can be improved) 
     708         zPr =   13.8_wp 
     709         zSc = 2432.0_wp 
     710 
     711         !! compute gamma mole 
     712         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 
     713         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 
     714 
     715         !! compute gamma 
     716         DO ji=2,jpi 
     717            DO jj=2,jpj 
    751718               ikt = mikt(ji,jj) 
    752719 
    753       !! Compute U and V at T points 
    754                zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) ) 
    755                zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) ) 
    756  
    757       !! compute ustar 
    758                zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 
    759                IF (zustar == 0._wp) THEN           ! only for kt = 1 I think 
    760                  gt = rn_gammat0 
    761                  gs = rn_gammas0 
     720               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think 
     721                  pgt = rn_gammat0 
     722                  pgs = rn_gammas0 
    762723               ELSE 
    763       !! compute Rc number (as done in zdfric.F90) 
    764                zcoef = 0.5 / fse3w(ji,jj,ikt) 
    765                !                                            ! shear of horizontal velocity 
    766                zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   & 
    767                   &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
    768                zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   & 
    769                   &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
    770                !                                            ! richardson number (minimum value set to zero) 
    771                zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 ) 
    772  
    773       !! compute Pr and Sc number (can be improved) 
    774                zPr =   13.8 
    775                zSc = 2432.0 
    776  
    777       !! compute gamma mole 
    778                zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0 
    779                zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0 
    780  
    781       !! compute bouyancy  
    782                zts(jp_tem) = ttbl(ji,jj) 
    783                zts(jp_sal) = stbl(ji,jj) 
    784                zdep        = fsdepw(ji,jj,ikt) 
    785                ! 
    786                CALL eos_rab( zts, zdep, zab ) 
     724                  !! compute Rc number (as done in zdfric.F90) 
     725                  zcoef = 0.5_wp / fse3w(ji,jj,ikt) 
     726                  !                                            ! shear of horizontal velocity 
     727                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  & 
     728                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
     729                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  & 
     730                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
     731                  !                                            ! richardson number (minimum value set to zero) 
     732                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 
     733 
     734                  !! compute bouyancy  
     735                  zts(jp_tem) = ttbl(ji,jj) 
     736                  zts(jp_sal) = stbl(ji,jj) 
     737                  zdep        = fsdepw(ji,jj,ikt) 
    787738                  ! 
    788       !! compute length scale  
    789                zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    790  
    791       !! compute Monin Obukov Length 
    792                ! Maximum boundary layer depth 
    793                zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001 
    794                ! Compute Monin obukhov length scale at the surface and Ekman depth: 
    795                zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln)) 
    796                zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 
    797  
    798       !! compute eta* (stability parameter) 
    799                zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0))) 
    800  
    801       !! compute the sublayer thickness 
    802                zhnu = 5 * znu / zustar 
    803       !! compute gamma turb 
    804                zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 
    805                &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn 
    806  
    807       !! compute gammats 
    808                gt = zustar / (zgturb + zgmolet) 
    809                gs = zustar / (zgturb + zgmoles) 
     739                  CALL eos_rab( zts, zdep, zab ) 
     740                  ! 
     741                  !! compute length scale  
     742                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     743 
     744                  !! compute Monin Obukov Length 
     745                  ! Maximum boundary layer depth 
     746                  zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp 
     747                  ! Compute Monin obukhov length scale at the surface and Ekman depth: 
     748                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 
     749                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 
     750 
     751                  !! compute eta* (stability parameter) 
     752                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 
     753 
     754                  !! compute the sublayer thickness 
     755                  zhnu = 5 * znu / zustar(ji,jj) 
     756 
     757                  !! compute gamma turb 
     758                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 
     759                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 
     760 
     761                  !! compute gammats 
     762                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 
     763                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 
    810764               END IF 
    811       END IF 
    812  
    813    END SUBROUTINE 
    814  
    815    SUBROUTINE sbc_isf_tbl( varin, varout, cptin ) 
     765            END DO 
     766         END DO 
     767         CALL lbc_lnk(pgt(:,:),'T',1.) 
     768         CALL lbc_lnk(pgs(:,:),'T',1.) 
     769      END SELECT 
     770      CALL wrk_dealloc( jpi,jpj, zustar ) 
     771      ! 
     772   END SUBROUTINE sbc_isf_gammats 
     773 
     774   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 
    816775      !!---------------------------------------------------------------------- 
    817776      !!                  ***  SUBROUTINE sbc_isf_tbl  *** 
    818777      !! 
    819       !! ** Purpose : compute mean T/S/U/V in the boundary layer  
    820       !! 
    821       !!---------------------------------------------------------------------- 
    822       REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin 
    823       REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout 
     778      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 
     779      !! 
     780      !!---------------------------------------------------------------------- 
     781      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin 
     782      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout 
     783      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out 
     784      ! 
     785      REAL(wp) :: ze3, zhk 
     786      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 
     787 
     788      INTEGER :: ji, jj, jk                  ! loop index 
     789      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl 
     790      !!---------------------------------------------------------------------- 
     791      ! allocation 
     792      CALL wrk_alloc( jpi,jpj, zhisf_tbl) 
    824793       
    825       CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out 
    826  
    827       REAL(wp) :: ze3, zhk 
    828       REAL(wp), DIMENSION(:,:), POINTER :: zikt 
    829  
    830       INTEGER :: ji,jj,jk 
    831       INTEGER :: ikt,ikb 
    832       INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb 
    833  
    834       CALL wrk_alloc( jpi,jpj, mkt, mkb  ) 
    835       CALL wrk_alloc( jpi,jpj, zikt ) 
    836  
    837       ! get first and last level of tbl 
    838       mkt(:,:) = misfkt(:,:) 
    839       mkb(:,:) = misfkb(:,:) 
    840  
    841       varout(:,:)=0._wp 
    842       DO jj = 2,jpj 
    843          DO ji = 2,jpi 
    844             IF (ssmask(ji,jj) == 1) THEN 
    845                ikt = mkt(ji,jj) 
    846                ikb = mkb(ji,jj) 
     794      ! initialisation 
     795      pvarout(:,:)=0._wp 
     796    
     797      SELECT CASE ( cd_ptin ) 
     798      CASE ( 'U' ) ! compute U in the top boundary layer at T- point  
     799         DO jj = 1,jpj 
     800            DO ji = 1,jpi 
     801               ikt = miku(ji,jj) ; ikb = miku(ji,jj) 
     802               ! thickness of boundary layer at least the top level thickness 
     803               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 
     804 
     805               ! determine the deepest level influenced by the boundary layer 
     806               DO jk = ikt+1, mbku(ji,jj) 
     807                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
     808               END DO 
     809               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     810 
     811               ! level fully include in the ice shelf boundary layer 
     812               DO jk = ikt, ikb - 1 
     813                  ze3 = fse3u_n(ji,jj,jk) 
     814                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     815               END DO 
     816 
     817               ! level partially include in ice shelf boundary layer  
     818               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     819               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     820            END DO 
     821         END DO 
     822         DO jj = 2,jpj 
     823            DO ji = 2,jpi 
     824               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 
     825            END DO 
     826         END DO 
     827         CALL lbc_lnk(pvarout,'T',-1.) 
     828       
     829      CASE ( 'V' ) ! compute V in the top boundary layer at T- point  
     830         DO jj = 1,jpj 
     831            DO ji = 1,jpi 
     832               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 
     833               ! thickness of boundary layer at least the top level thickness 
     834               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 
     835 
     836               ! determine the deepest level influenced by the boundary layer 
     837               DO jk = ikt+1, mbkv(ji,jj) 
     838                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
     839               END DO 
     840               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     841 
     842               ! level fully include in the ice shelf boundary layer 
     843               DO jk = ikt, ikb - 1 
     844                  ze3 = fse3v_n(ji,jj,jk) 
     845                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     846               END DO 
     847 
     848               ! level partially include in ice shelf boundary layer  
     849               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     850               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     851            END DO 
     852         END DO 
     853         DO jj = 2,jpj 
     854            DO ji = 2,jpi 
     855               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 
     856            END DO 
     857         END DO 
     858         CALL lbc_lnk(pvarout,'T',-1.) 
     859 
     860      CASE ( 'T' ) ! compute T in the top boundary layer at T- point  
     861         DO jj = 1,jpj 
     862            DO ji = 1,jpi 
     863               ikt = misfkt(ji,jj) 
     864               ikb = misfkb(ji,jj) 
    847865 
    848866               ! level fully include in the ice shelf boundary layer 
    849867               DO jk = ikt, ikb - 1 
    850868                  ze3 = fse3t_n(ji,jj,jk) 
    851                   IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    852                   IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) & 
    853                      &                                                       * r1_hisf_tbl(ji,jj) * ze3 
    854                   IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) & 
    855                      &                                                       * r1_hisf_tbl(ji,jj) * ze3 
     869                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    856870               END DO 
    857871 
    858872               ! level partially include in ice shelf boundary layer  
    859873               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
    860                IF (cptin == 'T') & 
    861                    &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
    862                IF (cptin == 'U') & 
    863                    &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 
    864                IF (cptin == 'V') & 
    865                    &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 
    866             END IF 
     874               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     875            END DO 
    867876         END DO 
    868       END DO 
    869  
    870       CALL wrk_dealloc( jpi,jpj, mkt, mkb )       
    871       CALL wrk_dealloc( jpi,jpj, zikt )  
    872  
    873       IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.) 
    874       IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.) 
    875  
     877      END SELECT 
     878 
     879      ! mask mean tbl value 
     880      pvarout(:,:) = pvarout(:,:) * ssisfmask(:,:) 
     881 
     882      ! deallocation 
     883      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )       
     884      ! 
    876885   END SUBROUTINE sbc_isf_tbl 
    877886       
    … …  
    889898      !! ** Action  :   phdivn   decreased by the runoff inflow 
    890899      !!---------------------------------------------------------------------- 
    891       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
    892       !! 
     900      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence 
     901      !  
    893902      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    894903      INTEGER  ::   ikt, ikb  
    895       INTEGER  ::   nk_isf 
    896       REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl 
    897       REAL(wp)     ::   zfact     ! local scalar 
     904      REAL(wp) ::   zhk 
     905      REAL(wp) ::   zfact     ! local scalar 
    898906      !!---------------------------------------------------------------------- 
    899907      ! 
    900908      zfact   = 0.5_wp 
    901909      ! 
    902       IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water 
     910      IF ( lk_vvl ) THEN     ! need to re compute level distribution of isf fresh water 
    903911         DO jj = 1,jpj 
    904912            DO ji = 1,jpi 
    … …  
    909917 
    910918               ! determine the deepest level influenced by the boundary layer 
    911                ! test on tmask useless ????? 
    912                DO jk = ikt, mbkt(ji,jj) 
    913                   IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
     919               DO jk = ikt+1, mbkt(ji,jj) 
     920                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) <  rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    914921               END DO 
    915922               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    … …  
    917924               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    918925 
    919                zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
     926               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
    920927               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
    921928            END DO 
    922929         END DO 
    923       END IF ! vvl case 
    924       ! 
     930      END IF  
     931      ! 
     932      !==   ice shelf melting distributed over several levels   ==! 
    925933      DO jj = 1,jpj 
    926934         DO ji = 1,jpi 
    … …  
    930938               DO jk = ikt, ikb - 1 
    931939                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 
    932                     &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
     940                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
    933941               END DO 
    934942               ! level partially include in ice shelf boundary layer  
    935943               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 
    936                   &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    937             !==   ice shelf melting mass distributed over several levels   ==! 
     944                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    938945         END DO 
    939946      END DO 
    940947      ! 
    941948   END SUBROUTINE sbc_isf_div 
    942                          
    943    FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti ) 
    944       !!---------------------------------------------------------------------- 
    945       !!                 ***  ROUTINE eos_init  *** 
    946       !! 
    947       !! ** Purpose :   Compute the in-situ temperature [Celcius] 
    948       !! 
    949       !! ** Method  :    
    950       !! 
    951       !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408 
    952       !!---------------------------------------------------------------------- 
    953       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius] 
    954       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu] 
    955       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar] 
    956       REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius] 
    957 !      REAL(wp) :: fsatg 
    958 !      REAL(wp) :: pfps, pfpt, pfphp  
    959       REAL(wp) :: zt, zs, zp, zh, zq, zxk 
    960       INTEGER  :: ji, jj            ! dummy loop indices 
    961       ! 
    962       CALL wrk_alloc( jpi,jpj, pti  ) 
    963       !  
    964       DO jj=1,jpj 
    965          DO ji=1,jpi 
    966             zh = ppress(ji,jj) 
    967 ! Theta1 
    968             zt = ptem(ji,jj) 
    969             zs = psal(ji,jj) 
    970             zp = 0.0 
    971             zxk= zh * fsatg( zs, zt, zp ) 
    972             zt = zt + 0.5 * zxk 
    973             zq = zxk 
    974 ! Theta2 
    975             zp = zp + 0.5 * zh 
    976             zxk= zh*fsatg( zs, zt, zp ) 
    977             zt = zt + 0.29289322 * ( zxk - zq ) 
    978             zq = 0.58578644 * zxk + 0.121320344 * zq 
    979 ! Theta3 
    980             zxk= zh * fsatg( zs, zt, zp ) 
    981             zt = zt + 1.707106781 * ( zxk - zq ) 
    982             zq = 3.414213562 * zxk - 4.121320344 * zq 
    983 ! Theta4 
    984             zp = zp + 0.5 * zh 
    985             zxk= zh * fsatg( zs, zt, zp ) 
    986             pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0 
    987          END DO 
    988       END DO 
    989       ! 
    990       CALL wrk_dealloc( jpi,jpj, pti  ) 
    991       ! 
    992    END FUNCTION tinsitu 
    993    ! 
    994    FUNCTION fsatg( pfps, pfpt, pfphp ) 
    995       !!---------------------------------------------------------------------- 
    996       !!                 ***  FUNCTION fsatg  *** 
    997       !! 
    998       !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1 
    999       !! 
    1000       !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408 
    1001       !!  
    1002       !! ** units      :   pressure        pfphp    decibars 
    1003       !!                   temperature     pfpt     deg celsius (ipts-68) 
    1004       !!                   salinity        pfps     (ipss-78) 
    1005       !!                   adiabatic       fsatg    deg. c/decibar 
    1006       !!---------------------------------------------------------------------- 
    1007       REAL(wp) :: pfps, pfpt, pfphp  
    1008       REAL(wp) :: fsatg 
    1009       ! 
    1010       fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         & 
    1011         &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    & 
    1012         &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             & 
    1013         &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         & 
    1014         &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5 
    1015       ! 
    1016     END FUNCTION fsatg 
    1017     !!====================================================================== 
     949   !!====================================================================== 
    1018950END MODULE sbcisf 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r9321 r14572  
    183183         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp 
    184184         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 
    185          rdivisf       = 0.0_wp 
    186185      END IF 
    187186      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r9321 r14572  
    132132            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg 
    133133               ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study) 
    134                rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp 
     134               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp 
    135135            END WHERE 
    136136         ELSE                                                        ! use SST as runoffs temperature 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r9321 r14572  
    120120      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices   
    121121      INTEGER  ::   ikt, ikb  
    122       INTEGER  ::   nk_isf 
    123122      REAL(wp) ::   zfact, z1_e3t, zdep 
    124       REAL(wp) ::   zalpha, zhk 
    125       REAL(wp) ::  zt_frz, zpress 
     123      REAL(wp) ::   zt_frz, zpress 
    126124      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    127125      !!---------------------------------------------------------------------- 
    … …  
    162160         ELSE                                         ! No restart or restart not found: Euler forward time stepping 
    163161            zfact = 1._wp 
    164             sbc_tsc(:,:,:) = 0._wp 
    165162            sbc_tsc_b(:,:,:) = 0._wp 
    166163         ENDIF 
    … …  
    226223      ! 
    227224      IF( nn_isf > 0 ) THEN 
    228          zfact = 0.5e0 
     225         zfact = 0.5_wp 
    229226         DO jj = 2, jpj 
    230227            DO ji = fs_2, fs_jpim1 
    … …  
    234231    
    235232               ! level fully include in the ice shelf boundary layer 
    236                ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst) 
    237233               ! sign - because fwf sign of evapo (rnf sign of precip) 
    238234               DO jk = ikt, ikb - 1 
    239                ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    240235               ! compute trend 
    241                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          & 
    242                      &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) 
    243                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          & 
    244                      &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 
     236                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                & 
     237                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
     238                     &           * r1_hisf_tbl(ji,jj) 
    245239               END DO 
    246240    
    247241               ! level partially include in ice shelf boundary layer  
    248                ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    249242               ! compute trend 
    250                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           & 
    251                   &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
    252                tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           & 
    253                   &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)  
     243               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 & 
     244                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
     245                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
     246 
    254247            END DO 
    255248         END DO 
Note: See TracChangeset for help on using the changeset viewer.