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 9855 – NEMO

Changeset 9855


Ignore:
Timestamp:
2018-06-28T16:33:17+02:00 (6 years ago)
Author:
mathiot
Message:

update sbcisf from trunk

Location:
branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/CONFIG/SHARED/field_def.xml

    r9163 r9855  
    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/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r6487 r9855  
    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/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r8046 r9855  
    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/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r8046 r9855  
    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                  !: net heat 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      REAL(wp)                     ::   zt_frz, zpress 
     98      CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file 
     99      CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file 
     100      CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale 
     101      INTEGER           ::   ios           ! Local integer output status for namelist read 
     102 
     103      REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 
     104      REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d 
     105      REAL(wp), DIMENSION(:,:  ), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)  
     106      ! 
     107      !!--------------------------------------------------------------------- 
     108      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, & 
     109                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
    100110      ! 
    101111      ! 
    … …  
    121131         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk 
    122132         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl 
    123          IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf  
    124133         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
    125134         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 
    131135         ! 
    132136         ! Allocate public variable 
    … …  
    185189             
    186190            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
    187             ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 
     191            ALLOCATE( sf_fwfisf(1), STAT=ierror ) 
    188192            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) ) 
    190193            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' ) 
    192194         END IF 
    193195          
    … …  
    206208 
    207209      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 
     210         ! allocation 
     211         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  ) 
    208212 
    209213         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 
    … …  
    229233 
    230234         ! compute salf and heat flux 
    231          IF (nn_isf == 1) THEN 
    232             ! realistic ice shelf formulation 
     235         SELECT CASE ( nn_isf ) 
     236         CASE ( 1 )    ! realistic ice shelf formulation 
    233237            ! compute T/S/U/V for the top boundary layer 
    234238            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
    235239            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') 
     240            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U') 
     241            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V') 
    238242            ! iom print 
    239243            CALL iom_put('ttbl',ttbl(:,:)) 
    … …  
    244248            CALL sbc_isf_cav (kt) 
    245249 
    246          ELSE IF (nn_isf == 2) THEN 
    247             ! Beckmann and Goosse parametrisation  
     250         CASE ( 2 )    ! Beckmann and Goosse parametrisation  
    248251            stbl(:,:)   = soce 
    249252            CALL sbc_isf_bg03(kt) 
    250253 
    251          ELSE IF (nn_isf == 3) THEN 
    252             ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
     254         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    253255            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)  
     256            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf flux from the isf (fwfisf <0 mean melting)  
    255257 
    256258            IF( lk_oasis) THEN 
    … …  
    294296            ENDIF 
    295297 
    296             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
     298            qisf(:,:)   = fwfisf(:,:) * rlfusisf              ! heat        flux 
    297299            stbl(:,:)   = soce 
    298300 
    299          ELSE IF (nn_isf == 4) THEN 
    300             ! specified fwf and heat flux forcing beneath the ice shelf 
     301         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf 
    301302            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 
     303            fwfisf(:,:) =   sf_fwfisf(1)%fnow(:,:,1)            ! fwf  flux from the isf (fwfisf <0 mean melting) 
    304304 
    305305            IF( lk_oasis) THEN 
    … …  
    343343            ENDIF 
    344344 
    345             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
    346             !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux 
     345            qisf(:,:)   = fwfisf(:,:) * rlfusisf              ! heat        flux 
    347346            stbl(:,:)   = soce 
    348347 
    349          END IF 
     348         END SELECT 
     349 
    350350         ! 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 ! 
     351         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 
     352         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 
     353         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 
     354         DO jj = 1,jpj 
     355            DO ji = 1,jpi 
     356               zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj)) 
     357            END DO 
     358         END DO 
     359         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 
    355360          
    356          ! salt effect already take into account in vertical advection 
    357          risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 
     361         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 
     362         risf_tsc(:,:,jp_sal) = 0.0_wp 
    358363 
    359364         ! output 
    … …  
    361366         IF( iom_use('fwfisf'  ) )   CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce ) 
    362367 
    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   
    366368         ! lbclnk 
    367369         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 
    368370         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 
    369          CALL lbc_lnk(fwfisf(:,:)   ,'T',1.) 
    370          CALL lbc_lnk(qisf(:,:)     ,'T',1.) 
     371         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.) 
     372         CALL lbc_lnk(qisf(:,:)           ,'T',1.) 
    371373 
    372374!============================================================================================================================================= 
    … …  
    405407!============================================================================================================================================= 
    406408 
    407          IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
     409         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    ! 
    408410            IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
    409411                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 
    … …  
    416418               risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 
    417419            END IF 
    418          ENDIF 
     420         END IF 
    419421         !  
     422         ! deallocation 
     423         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  ) 
    420424      END IF 
    421    
     425      !   
    422426  END SUBROUTINE sbc_isf 
     427 
    423428 
    424429  INTEGER FUNCTION sbc_isf_alloc() 
    … …  
    435440               &    STAT= sbc_isf_alloc ) 
    436441         ! 
    437          IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc ) 
     442         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc ) 
    438443         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 
    439444         ! 
    440       ENDIF 
     445      END IF 
    441446  END FUNCTION 
    442447 
    443448  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)  
     449      !!--------------------------------------------------------------------- 
     450      !!                  ***  ROUTINE sbc_isf_bg03  *** 
     451      !! 
     452      !! ** Purpose : add net heat and fresh water flux from ice shelf melting 
     453      !!          into the adjacent ocean 
     454      !! 
     455      !! ** Method  :   See reference 
     456      !! 
     457      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
     458      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170. 
     459      !!         (hereafter BG) 
     460      !! History : 
     461      !!         06-02  (C. Wang) Original code 
     462      !!---------------------------------------------------------------------- 
     463      INTEGER, INTENT ( in ) :: kt 
     464      ! 
     465      INTEGER  :: ji, jj, jk ! dummy loop index 
     466      INTEGER  :: ik         ! current level 
     467      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m 
     468      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m 
     469      REAL(wp) :: zt_frz     ! freezing point temperature at depth z 
     470      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth 
     471      !!---------------------------------------------------------------------- 
     472 
     473      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
     474      ! 
     475      DO ji = 1, jpi 
     476         DO jj = 1, jpj 
     477            ik = misfkt(ji,jj) 
     478            !! Initialize arrays to 0 (each step) 
     479            zt_sum = 0.e0_wp 
     480            IF ( ik > 1 ) THEN 
     481               ! 1. -----------the average temperature between 200m and 600m --------------------- 
     482               DO jk = misfkt(ji,jj),misfkb(ji,jj) 
     483                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 
     484                  ! after verif with UNESCO, wrong sign in BG eq. 2 
     485                  ! Calculate freezing temperature 
     486                  CALL eos_fzp(stbl(ji,jj), zt_frz, gdept_n(ji,jj,jk))  
     487                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
     488               END DO 
     489               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
     490               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 
     491               ! For those corresponding to zonal boundary     
     492               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
     493                           & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk) 
    501494              
    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') 
     495               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                   
     496               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
     497               !add to salinity trend 
     498            ELSE 
     499               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 
     500            END IF 
     501         END DO 
     502      END DO 
     503      ! 
     504      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03') 
     505      ! 
    512506  END SUBROUTINE sbc_isf_bg03 
    513507 
    … …  
    527521      INTEGER, INTENT(in)          ::   kt         ! ocean time step 
    528522      ! 
    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 
     523      INTEGER  ::   ji, jj     ! dummy loop indices 
     524      INTEGER  ::   nit 
    533525      REAL(wp) ::   zlamb1, zlamb2, zlamb3 
    534526      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 
     527      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zcfac 
     528      REAL(wp) ::   zeps = 1.e-20_wp         
     529      REAL(wp) ::   zerr 
     530      REAL(wp), DIMENSION(:,:), POINTER ::   ztfrz, zsfrz 
     531      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas  
     532      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b 
     533      LOGICAL  ::   lit 
    544534      !!--------------------------------------------------------------------- 
    545       ! 
    546       ! coeficient for linearisation of tfreez 
    547       zlamb1=-0.0575 
    548       zlamb2=0.0901 
    549       zlamb3=-7.61e-04 
     535      ! coeficient for linearisation of potential tfreez 
     536      ! Crude approximation for pressure (but commonly used) 
     537      IF ( ln_useCT ) THEN   ! linearisation from Jourdain et al. (2017) 
     538         zlamb1 =-0.0564_wp 
     539         zlamb2 = 0.0773_wp 
     540         zlamb3 =-7.8633e-8 * grav * rau0 
     541      ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015) 
     542         zlamb1 =-0.0573_wp 
     543         zlamb2 = 0.0832_wp 
     544         zlamb3 =-7.53e-8 * grav * rau0 
     545      ENDIF 
     546      ! 
    550547      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav') 
    551548      ! 
    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 
    567             ! 
    568          END DO 
     549      CALL wrk_alloc( jpi,jpj, ztfrz  , zsfrz  , zgammat, zgammas  ) 
     550      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 
     551 
     552      ! initialisation 
     553      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 
     554      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp     
     555      zfwflx (:,:) = 0.0_wp     ; zsfrz(:,:) = 0.0_wp 
     556 
     557      ! compute ice shelf melting 
     558      nit = 1 ; lit = .TRUE. 
     559      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine 
     560         SELECT CASE ( nn_isfblk ) 
     561         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 
     562            ! Calculate freezing temperature 
     563            CALL eos_fzp( stbl(:,:), ztfrz(:,:), risfdep(:,:) ) 
     564 
     565            ! compute gammat every where (2d) 
     566            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     567             
     568            ! compute upward heat flux zhtflx and upward water flux zwflx 
     569            DO jj = 1, jpj 
     570               DO ji = 1, jpi 
     571                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-ztfrz(ji,jj)) 
     572                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 
     573               END DO 
     574            END DO 
     575 
     576            ! Compute heat flux and upward fresh water flux 
     577            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     578            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     579 
     580         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
     581            ! compute gammat every where (2d) 
     582            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     583 
     584            ! compute upward heat flux zhtflx and upward water flux zwflx 
     585            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 
     586            DO jj = 1, jpj 
     587               DO ji = 1, jpi 
     588                  ! compute coeficient to solve the 2nd order equation 
     589                  zeps1 = rcp*rau0*zgammat(ji,jj) 
     590                  zeps2 = rlfusisf*rau0*zgammas(ji,jj) 
     591                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 
     592                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 
     593                  zeps6 = zeps4-ttbl(ji,jj) 
     594                  zeps7 = zeps4-tsurf 
     595                  zaqe  = zlamb1 * (zeps1 + zeps3) 
     596                  zaqer = 0.5_wp/MIN(zaqe,-zeps) 
     597                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2 
     598                  zcqe  = zeps2*stbl(ji,jj) 
     599                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe                
     600 
     601                  ! Presumably zdis can never be negative because gammas is very small compared to gammat 
     602                  ! compute s freeze 
     603                  zsfrz(ji,jj)=(-zbqe-SQRT(zdis))*zaqer 
     604                  IF ( zsfrz(ji,jj) < 0.0_wp ) zsfrz(ji,jj)=(-zbqe+SQRT(zdis))*zaqer 
     605 
     606                  ! compute t freeze (eq. 22) 
     607                  ztfrz(ji,jj)=zeps4+zlamb1*zsfrz(ji,jj) 
     608   
     609                  ! zfwflx is upward water flux 
     610                  ! zhtflx is upward heat flux (out of ocean) 
     611                  ! compute the upward water and heat flux (eq. 28 and eq. 29) 
     612                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz(ji,jj)-stbl(ji,jj)) / MAX(zsfrz(ji,jj),zeps) 
     613                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - ztfrz(ji,jj) )  
     614               END DO 
     615            END DO 
     616 
     617            ! compute heat and water flux 
     618            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     619            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     620 
     621         END SELECT 
     622 
     623         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 
     624         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE. 
     625         ELSE                            
     626            ! check total number of iteration 
     627            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
     628            ELSE                 ; nit = nit + 1 
     629            END IF 
     630 
     631            ! compute error between 2 iterations 
     632            ! if needed save gammat and compute zhtflx_b for next iteration 
     633            zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 
     634            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 
     635            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:) 
     636            END IF 
     637         END IF 
    569638      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 ) 
     639      ! 
     640      CALL iom_put('isftfrz'  , ztfrz * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
     641      CALL iom_put('isfsfrz'  , zsfrz * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
     642      CALL iom_put('isfgammat', zgammat) 
     643      CALL iom_put('isfgammas', zgammas) 
     644      !  
     645      CALL wrk_dealloc( jpi,jpj, ztfrz  , zgammat, zgammas  ) 
     646      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 
    682647      ! 
    683648      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav') 
    684  
     649      ! 
    685650   END SUBROUTINE sbc_isf_cav 
    686651 
    687    SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit ) 
     652   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 
    688653      !!---------------------------------------------------------------------- 
    689654      !! ** Purpose    : compute the coefficient echange for heat flux 
    … …  
    694659      !!                Jenkins et al., 2010, JPO, p2298-2312 
    695660      !!--------------------------------------------------------------------- 
    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 
     661      REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 
     662      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 
     663      ! 
     664      INTEGER  :: ikt                         
     665      INTEGER  :: ji, jj                     ! loop index 
     666      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity 
    702667      REAL(wp) :: zdku, zdkv                 ! U, V shear  
    703668      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number  
    … …  
    709674      REAL(wp) :: zcoef                      ! temporary coef 
    710675      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 
     676      REAL(wp) :: zeps = 1.0e-20_wp     
     677      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant 
     678      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 
    715679      REAL(wp), DIMENSION(2) :: zts, zab 
    716680      !!--------------------------------------------------------------------- 
    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) 
     681      CALL wrk_alloc( jpi,jpj, zustar ) 
     682      ! 
     683      SELECT CASE ( nn_gammablk ) 
     684      CASE ( 0 ) ! gamma is constant (specified in namelist) 
     685         !! ISOMIP formulation (Hunter et al, 2006) 
     686         pgt(:,:) = rn_gammat0 
     687         pgs(:,:) = rn_gammas0 
     688 
     689      CASE ( 1 ) ! gamma is assume to be proportional to u* 
     690         !! Jenkins et al., 2010, JPO, p2298-2312 
     691         !! Adopted by Asay-Davis et al. (2015) 
     692 
     693         !! compute ustar (eq. 24) 
     694         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
     695 
     696         !! Compute gammats 
     697         pgt(:,:) = zustar(:,:) * rn_gammat0 
     698         pgs(:,:) = zustar(:,:) * rn_gammas0 
     699       
     700      CASE ( 2 ) ! gamma depends of stability of boundary layer 
     701         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
     702         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
     703         !! compute ustar 
     704         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
     705 
     706         !! compute Pr and Sc number (can be improved) 
     707         zPr =   13.8_wp 
     708         zSc = 2432.0_wp 
     709 
     710         !! compute gamma mole 
     711         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 
     712         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 
     713 
     714         !! compute gamma 
     715         DO ji=2,jpi 
     716            DO jj=2,jpj 
    751717               ikt = mikt(ji,jj) 
    752718 
    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 
     719               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think 
     720                  pgt = rn_gammat0 
     721                  pgs = rn_gammas0 
    762722               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 ) 
     723                  !! compute Rc number (as done in zdfric.F90) 
     724                  zcoef = 0.5_wp / fse3w(ji,jj,ikt) 
     725                  !                                            ! shear of horizontal velocity 
     726                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  & 
     727                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
     728                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  & 
     729                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
     730                  !                                            ! richardson number (minimum value set to zero) 
     731                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 
     732 
     733                  !! compute bouyancy  
     734                  zts(jp_tem) = ttbl(ji,jj) 
     735                  zts(jp_sal) = stbl(ji,jj) 
     736                  zdep        = fsdepw(ji,jj,ikt) 
    787737                  ! 
    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) 
     738                  CALL eos_rab( zts, zdep, zab ) 
     739                  ! 
     740                  !! compute length scale  
     741                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     742 
     743                  !! compute Monin Obukov Length 
     744                  ! Maximum boundary layer depth 
     745                  zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp 
     746                  ! Compute Monin obukhov length scale at the surface and Ekman depth: 
     747                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 
     748                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 
     749 
     750                  !! compute eta* (stability parameter) 
     751                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 
     752 
     753                  !! compute the sublayer thickness 
     754                  zhnu = 5 * znu / zustar(ji,jj) 
     755 
     756                  !! compute gamma turb 
     757                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 
     758                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 
     759 
     760                  !! compute gammats 
     761                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 
     762                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 
    810763               END IF 
    811       END IF 
    812  
    813    END SUBROUTINE 
    814  
    815    SUBROUTINE sbc_isf_tbl( varin, varout, cptin ) 
     764            END DO 
     765         END DO 
     766         CALL lbc_lnk(pgt(:,:),'T',1.) 
     767         CALL lbc_lnk(pgs(:,:),'T',1.) 
     768      END SELECT 
     769      CALL wrk_dealloc( jpi,jpj, zustar ) 
     770      ! 
     771   END SUBROUTINE sbc_isf_gammats 
     772 
     773   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 
    816774      !!---------------------------------------------------------------------- 
    817775      !!                  ***  SUBROUTINE sbc_isf_tbl  *** 
    818776      !! 
    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 
     777      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 
     778      !! 
     779      !!---------------------------------------------------------------------- 
     780      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin 
     781      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout 
     782      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out 
     783      ! 
     784      REAL(wp) :: ze3, zhk 
     785      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 
     786 
     787      INTEGER :: ji, jj, jk                  ! loop index 
     788      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl 
     789      !!---------------------------------------------------------------------- 
     790      ! allocation 
     791      CALL wrk_alloc( jpi,jpj, zhisf_tbl) 
    824792       
    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) 
     793      ! initialisation 
     794      pvarout(:,:)=0._wp 
     795    
     796      SELECT CASE ( cd_ptin ) 
     797      CASE ( 'U' ) ! compute U in the top boundary layer at T- point  
     798         DO jj = 1,jpj 
     799            DO ji = 1,jpi 
     800               ikt = miku(ji,jj) ; ikb = miku(ji,jj) 
     801               ! thickness of boundary layer at least the top level thickness 
     802               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 
     803 
     804               ! determine the deepest level influenced by the boundary layer 
     805               DO jk = ikt+1, mbku(ji,jj) 
     806                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
     807               END DO 
     808               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     809 
     810               ! level fully include in the ice shelf boundary layer 
     811               DO jk = ikt, ikb - 1 
     812                  ze3 = fse3u_n(ji,jj,jk) 
     813                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     814               END DO 
     815 
     816               ! level partially include in ice shelf boundary layer  
     817               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     818               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     819            END DO 
     820         END DO 
     821         DO jj = 2,jpj 
     822            DO ji = 2,jpi 
     823               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 
     824            END DO 
     825         END DO 
     826         CALL lbc_lnk(pvarout,'T',-1.) 
     827       
     828      CASE ( 'V' ) ! compute V in the top boundary layer at T- point  
     829         DO jj = 1,jpj 
     830            DO ji = 1,jpi 
     831               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 
     832               ! thickness of boundary layer at least the top level thickness 
     833               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 
     834 
     835               ! determine the deepest level influenced by the boundary layer 
     836               DO jk = ikt+1, mbkv(ji,jj) 
     837                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
     838               END DO 
     839               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     840 
     841               ! level fully include in the ice shelf boundary layer 
     842               DO jk = ikt, ikb - 1 
     843                  ze3 = fse3v_n(ji,jj,jk) 
     844                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     845               END DO 
     846 
     847               ! level partially include in ice shelf boundary layer  
     848               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     849               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     850            END DO 
     851         END DO 
     852         DO jj = 2,jpj 
     853            DO ji = 2,jpi 
     854               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 
     855            END DO 
     856         END DO 
     857         CALL lbc_lnk(pvarout,'T',-1.) 
     858 
     859      CASE ( 'T' ) ! compute T in the top boundary layer at T- point  
     860         DO jj = 1,jpj 
     861            DO ji = 1,jpi 
     862               ikt = misfkt(ji,jj) 
     863               ikb = misfkb(ji,jj) 
    847864 
    848865               ! level fully include in the ice shelf boundary layer 
    849866               DO jk = ikt, ikb - 1 
    850867                  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 
     868                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    856869               END DO 
    857870 
    858871               ! level partially include in ice shelf boundary layer  
    859872               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 
    867          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  
     873               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     874            END DO 
     875         END DO 
     876      END SELECT 
     877 
     878      ! mask mean tbl value 
     879      pvarout(:,:) = pvarout(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     880 
     881      ! deallocation 
     882      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )       
     883      ! 
    876884   END SUBROUTINE sbc_isf_tbl 
    877885       
    … …  
    889897      !! ** Action  :   phdivn   decreased by the runoff inflow 
    890898      !!---------------------------------------------------------------------- 
    891       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
    892       !! 
     899      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence 
     900      !  
    893901      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    894902      INTEGER  ::   ikt, ikb  
    895       INTEGER  ::   nk_isf 
    896       REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl 
    897       REAL(wp)     ::   zfact     ! local scalar 
     903      REAL(wp) ::   zhk 
     904      REAL(wp) ::   zfact     ! local scalar 
    898905      !!---------------------------------------------------------------------- 
    899906      ! 
    900907      zfact   = 0.5_wp 
    901908      ! 
    902       IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water 
     909      IF ( lk_vvl ) THEN     ! need to re compute level distribution of isf fresh water 
    903910         DO jj = 1,jpj 
    904911            DO ji = 1,jpi 
    … …  
    909916 
    910917               ! 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 
     918               DO jk = ikt+1, mbkt(ji,jj) 
     919                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) <  rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    914920               END DO 
    915921               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    … …  
    917923               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    918924 
    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 
     925               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
    920926               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
    921927            END DO 
    922928         END DO 
    923       END IF ! vvl case 
    924       ! 
     929      END IF  
     930      ! 
     931      !==   ice shelf melting distributed over several levels   ==! 
    925932      DO jj = 1,jpj 
    926933         DO ji = 1,jpi 
    … …  
    930937               DO jk = ikt, ikb - 1 
    931938                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 
    932                     &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
     939                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
    933940               END DO 
    934941               ! level partially include in ice shelf boundary layer  
    935942               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   ==! 
     943                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    938944         END DO 
    939945      END DO 
    940946      ! 
    941947   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     !!====================================================================== 
     948   !!====================================================================== 
    1018949END MODULE sbcisf 

  • branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r9321 r9855  
    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/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r9321 r9855  
    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/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r9321 r9855  
    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.