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 12077 for NEMO/branches – NEMO

Changeset 12077 for NEMO/branches


Ignore:
Timestamp:
2019-12-05T18:41:39+01:00 (4 years ago)
Author:
mathiot
Message:

include ENHANCE-02_ISF_nemo in UKMO merge branch

Location:
NEMO/branches/2019/UKMO_MERGE_2019
Files:
2 deleted
48 edited
2 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml

    r12068 r12077  
    264264 
    265265          <!-- * variable related to ice shelf forcing * --> 
    266           <field id="isftfrz_par"     long_name="fzp temperature at ocean/isf interface"        unit="degC"     /> 
    267           <field id="isftfrz_cav"     long_name="fzp temperature at ocean/isf interface"        unit="degC"     /> 
    268           <field id="fwfisf_cav"      long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
    269           <field id="fwfisf_par"      long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
     266          <field id="isftfrz_cav"     long_name="freezing point temperature at ocean/isf interface"                unit="degC"     /> 
     267          <field id="isftfrz_par"     long_name="freezing point temperature in the parametrization boundary layer" unit="degC"     /> 
     268          <field id="fwfisf_cav"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  /> 
     269          <field id="fwfisf_par"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  /> 
    270270          <field id="qoceisf_cav"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     /> 
    271271          <field id="qoceisf_par"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     /> 
    272272          <field id="qlatisf_cav"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     /> 
    273273          <field id="qlatisf_par"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     /> 
    274           <field id="qhcisf_cav"      long_name="Ice shelf heat content flux"                  unit="W/m2"     /> 
    275           <field id="qhcisf_par"      long_name="Ice shelf heat content flux"                  unit="W/m2"     /> 
    276           <field id="fwfisf3d_cav"    long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" /> 
    277           <field id="fwfisf3d_par"    long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" /> 
     274          <field id="qhcisf_cav"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     /> 
     275          <field id="qhcisf_par"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     /> 
     276          <field id="fwfisf3d_cav"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" /> 
     277          <field id="fwfisf3d_par"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" /> 
    278278          <field id="qoceisf3d_cav"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
    279279          <field id="qoceisf3d_par"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
    280280          <field id="qlatisf3d_cav"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
    281281          <field id="qlatisf3d_par"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
    282           <field id="qhcisf3d_cav"    long_name="Ice shelf heat content flux"                  unit="W/m2"     grid_ref="grid_T_3D" /> 
    283           <field id="qhcisf3d_par"    long_name="Ice shelf heat content flux"                  unit="W/m2"     grid_ref="grid_T_3D" /> 
    284           <field id="ttbl_cav"        long_name="temperature in the tracer sample depth      "  unit="C"        /> 
    285           <field id="ttbl_par"        long_name="temperature in the tracer sample depth      "  unit="C"        /> 
    286           <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting        "  unit="C"        /> 
    287           <field id="isfthermald_par" long_name="thermal driving of ice shelf melting        "  unit="C"        /> 
    288           <field id="isfgammat"       long_name="transfert coefficient for isf (temperature) "  unit="m/s"      /> 
    289           <field id="isfgammas"       long_name="transfert coefficient for isf (salinity)    "  unit="m/s"      /> 
    290           <field id="stbl"            long_name="salinity in the Losh tbl                    "  unit="PSU"      /> 
    291           <field id="utbl"            long_name="zonal current in the Losh tbl at T point    "  unit="m/s"      /> 
    292           <field id="vtbl"            long_name="merid current in the Losh tbl at T point    "  unit="m/s"      /> 
    293           <field id="isfustar"        long_name="ustar at T point used in ice shelf melting  "  unit="m/s"      /> 
     282          <field id="qhcisf3d_cav"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" /> 
     283          <field id="qhcisf3d_par"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" /> 
     284          <field id="ttbl_cav"        long_name="temperature in Losch tbl"                      unit="degC"     /> 
     285          <field id="ttbl_par"        long_name="temperature in the parametrisation boundary layer" unit="degC" /> 
     286          <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting"          unit="degC"     /> 
     287          <field id="isfthermald_par" long_name="thermal driving of ice shelf melting"          unit="degC"     /> 
     288          <field id="isfgammat"       long_name="Ice shelf heat-transfert velocity"             unit="m/s"      /> 
     289          <field id="isfgammas"       long_name="Ice shelf salt-transfert velocity"             unit="m/s"      /> 
     290          <field id="stbl"            long_name="salinity in the Losh tbl"                      unit="1e-3"     /> 
     291          <field id="utbl"            long_name="zonal current in the Losh tbl at T point"      unit="m/s"      /> 
     292          <field id="vtbl"            long_name="merid current in the Losh tbl at T point"      unit="m/s"      /> 
     293          <field id="isfustar"        long_name="ustar at T point used in ice shelf melting  unit="m/s"      /> 
    294294          <field id="qconisf"         long_name="Conductive heat flux through the ice shelf"    unit="W/m2"     /> 
    295295 
  • NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/namelist_ref

    r12068 r12077  
    439439   ! ---------------- ice shelf load ------------------------------- 
    440440   ! 
    441    cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc) 
     441   cn_isfload = 'uniform'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc) 
     442      rn_isfload_T = -1.9 
     443      rn_isfload_S = 34.4 
    442444   ! 
    443445   ! ---------------- ice shelf melt formulation ------------------------------- 
     
    452454         cn_isfcav_mlt = '3eq'   ! ice shelf melting formulation (spe/2eq/3eq/oasis) 
    453455         !                       ! spe = fwfisf is read from a forcing field 
    454          !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006) 
    455          !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 
     456         !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006 for a short description) 
     457         !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2016 for a short description) 
    456458         !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf 
    457459         !              !  cn_isfcav_mlt = 2eq or 3eq cases: 
    458          cn_gammablk = 'ad15'    ! scheme to compute gammat/s (spe,ad15,hj99) 
    459          !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
    460          !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999) 
    461          rn_gammat0  = 1.4e-2    ! gammat coefficient used in blk formula 
    462          rn_gammas0  = 4.e-4    ! gammas coefficient used in blk formula 
     460         cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99) 
     461         !                       ! spe      = constant transfert velocity (rn_gammat0, rn_gammas0) 
     462         !                       ! vel      = velocity dependent transfert velocity (u* * gammat/s) (Asay-Davis et al. 2016 for a short description) 
     463         !                       ! vel_stab = velocity and stability dependent transfert coeficient (Holland et al. 1999 for a complete description) 
     464         rn_gammat0  = 1.4e-2    ! gammat coefficient used in spe, vel and vel_stab gamma computation method 
     465         rn_gammas0  = 4.0e-4    ! gammas coefficient used in spe, vel and vel_stab gamma computation method 
    463466         ! 
    464467         rn_htbl     =  30.      ! thickness of the top boundary layer    (Losh et al. 2008) 
  • NEMO/branches/2019/UKMO_MERGE_2019/cfgs/WED025/EXPREF/namelist_cfg

    r11844 r12077  
    211211!----------------------------------------------------------------------- 
    212212   ! 
    213    ! ---------------- ice shelf load ------------------------------- 
    214    ! 
    215    cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc) 
    216    ! 
    217213   ! ---------------- ice shelf melt formulation ------------------------------- 
    218214   ! 
     
    229225         !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf 
    230226         !              !  cn_isfcav_mlt = 2eq or 3eq cases: 
    231          cn_gammablk = 'ad15'    ! scheme to compute gammat/s (spe,ad15,hj99) 
     227         cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99) 
    232228         !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
    233229         !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999) 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/BDY/bdyvol.F90

    r12068 r12077  
    1414   USE bdy_oce        ! ocean open boundary conditions 
    1515   USE sbc_oce        ! ocean surface boundary conditions 
     16   USE isf_oce, ONLY : fwfisf_cav, fwfisf_par  ! ice shelf 
    1617   USE dom_oce        ! ocean space and time domain  
    1718   USE phycst         ! physical constants 
    18    USE isf            ! ice shelf 
    1919   ! 
    2020   USE in_out_manager ! I/O manager 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diahsb.F90

    r12068 r12077  
    1717   USE phycst         ! physical constants 
    1818   USE sbc_oce        ! surface thermohaline fluxes 
     19   USE isf_oce        ! ice shelf fluxes 
    1920   USE sbcrnf         ! river runoff 
    20    USE isf            ! ice shelves 
    2121   USE domvvl         ! vertical scale factors 
    2222   USE traqsr         ! penetrative solar radiation 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diawri.F90

    r12068 r12077  
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce            ! ocean dynamics and tracers  
    28    USE isf 
     28   USE isf_oce 
    2929   USE isfcpl 
    3030   USE dom_oce        ! ocean space and time domain 
     
    911911      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity 
    912912      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height 
     913 
    913914      IF ( ln_isf ) THEN 
    914915         IF (ln_isfcav_mlt) THEN 
     
    918919            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity 
    919920            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity 
     921            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 ) 
    920922         END IF 
    921923         IF (ln_isfpar_mlt) THEN 
     
    926928            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity 
    927929            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity 
     930            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 ) 
    928931         END IF 
    929932      END IF 
    930933 
    931       IF ( ln_isf ) THEN 
    932          IF (ln_isfcav_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 ) 
    933          IF (ln_isfpar_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 ) 
    934       END IF 
    935        
    936934      IF( ALLOCATED(ahtu) ) THEN 
    937935         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DOM/domwri.F90

    r12068 r12077  
    1616   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 
    1717   !!---------------------------------------------------------------------- 
    18    USE isf             ! ice shelf 
     18   ! 
    1919   USE dom_oce         ! ocean space and time domain 
    2020   USE phycst ,   ONLY :   rsmall 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/divhor.F90

    r12068 r12077  
    1919   !!---------------------------------------------------------------------- 
    2020   USE oce             ! ocean dynamics and tracers 
    21    USE isf 
    22    USE isfutils 
    2321   USE dom_oce         ! ocean space and time domain 
    24    USE sbc_oce, ONLY : ln_rnf ! surface boundary condition: ocean 
    25    USE sbcrnf          ! river runoff  
    26    USE isfhdiv         ! ice shelf 
     22   USE sbc_oce, ONLY : ln_rnf      ! river runoff 
     23   USE sbcrnf , ONLY : sbc_rnf_div ! river runoff  
     24   USE isf_oce, ONLY : ln_isf      ! ice shelf 
     25   USE isfhdiv, ONLY : isf_hdiv    ! ice shelf 
    2726#if defined key_asminc    
    2827   USE asminc          ! Assimilation increment 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynatf.F90

    r12068 r12077  
    4141   USE trddyn         ! trend manager: dynamics 
    4242   USE trdken         ! trend manager: kinetic energy 
    43    USE isf       , ONLY: ln_isf     ! ice shelf 
    44    USE isfdynatf , ONLY: isf_dynatf ! ice shelf  
     43   USE isf_oce   , ONLY: ln_isf     ! ice shelf 
     44   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine  
    4545   ! 
    4646   USE in_out_manager ! I/O manager 
     
    223223            ! PM: we could probably define a generic subroutine to do the in depth correction 
    224224            !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...) 
     225            !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf) 
    225226            IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, atfp * rdt ) 
    226227            ! 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynhpg.F90

    r12068 r12077  
    3131   !!---------------------------------------------------------------------- 
    3232   USE oce             ! ocean dynamics and tracers 
    33    USE isf             ! ice shelf  (risfload variable) 
    34    USE isfload         ! ice shelf  (isf_load routine ) 
     33   USE isf_oce , ONLY : risfload  ! ice shelf  (risfload variable) 
     34   USE isfload , ONLY : isf_load  ! ice shelf  (isf_load routine ) 
    3535   USE sbc_oce         ! surface variable (only for the flag with ice shelf) 
    3636   USE dom_oce         ! ocean space and time domain 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynspg_ts.F90

    r12068 r12077  
    3131   USE dom_oce         ! ocean space and time domain 
    3232   USE sbc_oce         ! surface boundary condition: ocean 
     33   USE isf_oce         ! ice shelf variable (fwfisf) 
    3334   USE zdf_oce         ! vertical physics: variables 
    3435   USE zdfdrg          ! vertical physics: top/bottom drag coef. 
    35    USE isf             ! ice shelf variable (fwfisf) 
    36    USE isfutils 
    3736   USE sbcapr          ! surface boundary condition: atmospheric pressure 
    3837   USE dynadv    , ONLY: ln_dynadv_vec 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/sshwzv.F90

    r12068 r12077  
    1919   !!---------------------------------------------------------------------- 
    2020   USE oce            ! ocean dynamics and tracers variables 
    21    USE isf            ! ice shelf 
    22    USE isfutils 
     21   USE isf_oce        ! ice shelf 
    2322   USE dom_oce        ! ocean space and time domain variables  
    2423   USE sbc_oce        ! surface boundary condition: ocean 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcav.F90

    r12068 r12077  
    1313   !!   isf_cav       : update ice shelf melting under ice shelf 
    1414   !!---------------------------------------------------------------------- 
    15    USE isf            ! ice shelf public variables 
     15   USE isf_oce        ! ice shelf public variables 
    1616   ! 
    1717   USE isfrst   , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine 
     
    2525   USE par_oce  , ONLY: jpi,jpj                ! ocean space and time domain 
    2626   USE phycst   , ONLY: grav,rau0,r1_rau0_rcp  ! physical constants 
    27    USE eosbn2   , ONLY: l_useCT                ! l_useCT 
     27   USE eosbn2   , ONLY: ln_teos10               ! use ln_teos10 or not 
    2828   ! 
    2929   USE in_out_manager ! I/O manager 
     
    8585      ! 
    8686      ! initialisation 
    87       IF (TRIM(cn_gammablk) == 'HJ99' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rau0_rcp ! last time step total heat fluxes (to speed up convergence) 
     87      IF (TRIM(cn_gammablk) == 'vel_stab' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rau0_rcp ! last time step total heat fluxes (to speed up convergence) 
    8888      ! 
    8989      ! compute ice shelf melting 
     
    104104         ! define if we need to iterate 
    105105         SELECT CASE ( cn_gammablk ) 
    106          CASE ( 'spe','ad15' ) 
     106         CASE ( 'spe','vel' ) 
    107107            ! no convergence needed 
    108108            lit = .FALSE. 
    109          CASE ( 'hj99' ) 
     109         CASE ( 'vel_stab' ) 
    110110            ! compute error between 2 iterations 
    111111            zerr = MAXVAL(ABS(zqoce(:,:) - zqoce_b(:,:))) 
     
    113113            ! define if iteration needed 
    114114            IF (nit >= 100) THEN              ! too much iteration 
    115                CALL ctl_stop( 'STOP', 'isf_cav: HJ99 gamma formulation had too many iterations ...' ) 
     115               CALL ctl_stop( 'STOP', 'isf_cav: vel_stab gamma formulation had too many iterations ...' ) 
    116116            ELSE IF ( zerr <= 0.01_wp ) THEN  ! convergence is achieve 
    117117               lit = .FALSE. 
     
    212212         ! coeficient for linearisation of potential tfreez 
    213213         ! Crude approximation for pressure (but commonly used) 
    214          IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017) 
     214         IF ( ln_teos10 ) THEN   ! linearisation from Jourdain et al. (2017) 
    215215            risf_lamb1 =-0.0564_wp 
    216216            risf_lamb2 = 0.0773_wp 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavgam.F90

    r12068 r12077  
    1010   !!   isfcav_gammats       : compute exchange coeficient gamma  
    1111   !!---------------------------------------------------------------------- 
    12    USE isf 
     12   USE isf_oce 
    1313   USE isfutils, ONLY: debug 
    1414   USE isftbl  , ONLY: isf_tbl 
     
    4646      !! 
    4747      !! ** Method     : select the gamma formulation 
    48       !!                 3 method available (cst, AD15 and HJ99) 
     48      !!                 3 method available (cst, vel and vel_stab) 
    4949      !!--------------------------------------------------------------------- 
    5050      !!-------------------------- OUT ------------------------------------- 
     
    6666         ! gamma is constant (specified in namelist) 
    6767         ! nothing to do 
    68       CASE ('ad15', 'hj99') 
     68      CASE ('vel', 'vel_stab') 
    6969         ! compute velocity in tbl 
    7070         CALL isf_tbl(Kmm, uu(:,:,:,Kmm) ,zutbl(:,:),'U', miku, rhisf_tbl_cav) 
     
    9090         pgt(:,:) = rn_gammat0 
    9191         pgs(:,:) = rn_gammas0 
    92       CASE ( 'ad15' ) ! gamma is proportional to u* 
    93          CALL gammats_AD15 (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs ) 
    94       CASE ( 'hj99' ) ! gamma depends of stability of boundary layer and u* 
    95          CALL gammats_HJ99 (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs ) 
     92      CASE ( 'vel' ) ! gamma is proportional to u* 
     93         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs ) 
     94      CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u* 
     95         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs ) 
    9696      CASE DEFAULT 
    9797         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)') 
     
    116116   !!----------------------------------------------------------------------------------------------------- 
    117117   ! 
    118    SUBROUTINE gammats_AD15(putbl, pvtbl, pCd, pke2, &   ! <<== in 
     118   SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, &   ! <<== in 
    119119      &                                  pgt, pgs   )   ! ==>> out gammats [m/s] 
    120120      !!---------------------------------------------------------------------- 
     
    123123      !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar ) 
    124124      !! 
    125       !! ** Reference  : Jenkins et al., 2010, JPO, p2298-2312 
    126       !!                 Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016 
     125      !! ** Reference  : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016 
    127126      !!--------------------------------------------------------------------- 
    128127      !!-------------------------- OUT ------------------------------------- 
     
    146145      CALL iom_put('isfustar',zustar(:,:)) 
    147146      ! 
    148    END SUBROUTINE gammats_AD15 
    149  
    150    SUBROUTINE gammats_HJ99( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in 
    151       &                                                                 pgt  , pgs    )  ! ==>> out gammats [m/s] 
     147   END SUBROUTINE gammats_vel 
     148 
     149   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in 
     150      &                                                                     pgt  , pgs    )  ! ==>> out gammats [m/s] 
    152151      !!---------------------------------------------------------------------- 
    153152      !! ** Purpose    : compute the coefficient echange coefficient  
     
    186185      ! 
    187186      ! compute ustar 
    188       zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + r_ke0_top ) ) 
     187      zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) ) 
    189188      ! 
    190189      ! output ustar 
     
    254253      END DO 
    255254 
    256    END SUBROUTINE gammats_HJ99 
     255   END SUBROUTINE gammats_vel_stab 
    257256 
    258257END MODULE isfcavgam 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavmlt.F90

    r12068 r12077  
    1212   !!---------------------------------------------------------------------- 
    1313 
    14    USE isf                      ! ice shelf 
     14   USE isf_oce                  ! ice shelf 
    1515   USE isftbl , ONLY: isf_tbl   ! ice shelf depth average 
    1616   USE isfutils,ONLY: debug     ! debug subroutine 
     
    7777      ! 
    7878      IF (ln_isfdebug) THEN 
    79          CALL debug( 'isfcav_mlt:', pqhc (:,:) ) 
    80          CALL debug( 'isfcav_mlt:', pqoce(:,:) ) 
    81          CALL debug( 'isfcav_mlt:', pqfwf(:,:) ) 
     79         CALL debug( 'isfcav_mlt qhc  :', pqhc (:,:) ) 
     80         CALL debug( 'isfcav_mlt qoce :', pqoce(:,:) ) 
     81         CALL debug( 'isfcav_mlt qfwf :', pqfwf(:,:) ) 
    8282      END IF 
    8383      ! 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcpl.F90

    r12068 r12077  
    1212   !!   isfrst : read/write iceshelf variables in/from restart 
    1313   !!---------------------------------------------------------------------- 
    14    USE isf                              ! ice shelf variable 
     14   USE isf_oce                          ! ice shelf variable 
    1515   USE isfutils, ONLY : debug 
    1616   USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine 
     
    2525   PRIVATE 
    2626 
    27    PUBLIC isfcpl_rst_write, isfcpl_init  ! iceshelf restart read and write  
    28    PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons 
     27   PUBLIC isfcpl_rst_write, isfcpl_init                    ! iceshelf restart read and write  
     28   PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons  ! iceshelf correction for ssh, tra, dyn and conservation  
    2929 
    3030   TYPE isfcons 
     
    4848   SUBROUTINE isfcpl_init(Kbb, Kmm, Kaa) 
    4949      !!--------------------------------------------------------------------- 
     50      !!                   ***  ROUTINE iscpl_init  *** 
     51      !!  
     52      !! ** Purpose : correct ocean state for new wet cell and horizontal divergence  
     53      !!              correction for the dynamical adjustement 
     54      !! 
     55      !! ** Action : - compute ssh on new wet cell 
     56      !!             - compute T/S on new wet cell 
     57      !!             - compute horizontal divergence correction as a volume flux 
     58      !!             - compute the T/S/vol correction increment to keep trend to 0 
    5059      !! 
    5160      !!--------------------------------------------------------------------- 
     
    480489      INTEGER  ::   jisf                              ! start, end and current position in the increment array 
    481490      INTEGER  ::   ingb, ifind                       ! 0/1 target found or need to be found  
    482       INTEGER  ::   ingb, ifind, nisfl_area           ! global number of cell concerned by the wet->dry case  
     491      INTEGER  ::   nisfl_area                        ! global number of cell concerned by the wet->dry case  
    483492      INTEGER, DIMENSION(jpnij) :: nisfl              ! local  number of cell concerned by the wet->dry case 
    484493      ! 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdiags.F90

    r12068 r12077  
    11MODULE isfdiags 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  sbcisf  *** 
    4    !! Surface module :  update surface ocean boundary condition under ice 
    5    !!                   shelf 
     3   !!                       ***  MODULE  isfdiags  *** 
     4   !! ice shelf diagnostics module :  manage the 2d and 3d flux outputs from the ice shelf module 
    65   !!====================================================================== 
    76   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav 
     
    1615   USE in_out_manager ! I/O manager 
    1716   USE dom_oce 
    18    USE isf            ! ice shelf variable 
    19    USE isfutils 
     17   USE isf_oce        ! ice shelf variable 
    2018   USE iom            !  
    2119 
     
    3836      !!                  ***  ROUTINE isf_diags_flx *** 
    3937      !! 
    40       !! ** Purpose :  
     38      !! ** Purpose : manage the 2d and 3d flux outputs of the ice shelf module 
     39      !!              fwf, latent heat flux, heat content flux, oce->ice heat flux 
    4140      !! 
    4241      !!---------------------------------------------------------------------- 
     
    5958      cvarqhc  = 'qhcisf_'//cdisf  ; cvarqhc3d  = 'qhcisf3d_'//cdisf 
    6059      ! 
    61       ! output 2d melt rate, latent heat and heat content flux 
    62       CALL iom_put( TRIM(cvarqfwf), pqfwf(:,:) )   ! isf mass flux (opposite sign) 
    63       CALL iom_put( TRIM(cvarqoce), pqoce(:,:) )   ! isf oce to ice   flux  (cpo*gt*dt) 
    64       CALL iom_put( TRIM(cvarqlat), pqlat(:,:) )   ! isf oce to ice   flux  (cpo*gt*dt) 
    65       CALL iom_put( TRIM(cvarqhc) , pqhc (:,:) )   ! isf heat content flux  (cpo*fwf*Tfrz) 
     60      ! output 2d melt rate, latent heat and heat content flux from the injected water 
     61      CALL iom_put( TRIM(cvarqfwf), pqfwf(:,:) )   ! mass         flux ( >0 out ) 
     62      CALL iom_put( TRIM(cvarqoce), pqoce(:,:) )   ! oce to ice   flux ( >0 out ) 
     63      CALL iom_put( TRIM(cvarqlat), pqlat(:,:) )   ! latent heat  flux ( >0 out ) 
     64      CALL iom_put( TRIM(cvarqhc) , pqhc (:,:) )   ! heat content flux ( >0 out ) 
    6665      ! 
    6766      ! output 3d Diagnostics 
     
    7776      !!                  ***  ROUTINE isf_diags_2dto3d *** 
    7877      !! 
    79       !! ** Purpose :  
     78      !! ** Purpose : compute the 3d flux outputs as they are injected into NEMO  
     79      !!              (ie uniformaly spread into the top boundary layer or parametrisation layer) 
    8080      !! 
    8181      !!---------------------------------------------------------------------- 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdynatf.F90

    r12068 r12077  
    1111   !!------------------------------------------------------------------------- 
    1212 
    13    USE isf 
     13   USE isf_oce 
    1414 
    1515   USE phycst , ONLY: r1_rau0         ! physical constant 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfhdiv.F90

    r12068 r12077  
    11MODULE isfhdiv 
     2   !!====================================================================== 
     3   !!                       ***  MODULE  isfhdiv  *** 
     4   !! ice shelf horizontal divergence module :  update the horizontal divergence 
     5   !!                   with the ice shelf melt and coupling correction 
     6   !!====================================================================== 
     7   !! History :  4.0  !  2019-09  (P. Mathiot) Original code 
     8   !!---------------------------------------------------------------------- 
    29 
    3    USE isf                    ! ice shelf 
    4    USE isfutils  
     10   !!---------------------------------------------------------------------- 
     11   !!   isf_hdiv    : update the horizontal divergence with the ice shelf  
     12   !!                 melt and coupling correction 
     13   !!---------------------------------------------------------------------- 
     14 
     15   USE isf_oce                ! ice shelf 
     16 
    517   USE dom_oce                ! time and space domain 
    618   USE phycst , ONLY: r1_rau0 ! physical constant 
     
    4052         IF ( ln_isfcpl .AND. kt /= 0 ) THEN 
    4153            ! 
    42             ! correct divergence only for the first time step 
     54            ! Dynamical stability at start up after change in under ice shelf cavity geometry is achieve by correcting the divergence. 
     55            ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping  
     56            ! the same as at the end of the latest time step. So correction need to be apply at nit000 (euler time step) and 
     57            ! half of it at nit000+1 (leap frog time step). 
    4358            IF ( kt == nit000   ) CALL isf_hdiv_cpl(Kmm, risfcpl_vol       , phdiv) 
    4459            IF ( kt == nit000+1 ) CALL isf_hdiv_cpl(Kmm, risfcpl_vol*0.5_wp, phdiv) 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfload.F90

    r12068 r12077  
    1111   !!---------------------------------------------------------------------- 
    1212 
    13    USE isf, ONLY: cn_isfload  ! ice shelf variables 
     13   USE isf_oce, ONLY: cn_isfload, rn_isfload_T, rn_isfload_S ! ice shelf variables 
    1414 
    1515   USE dom_oce, ONLY: e3w, gdept, risfdep, mikt     ! vertical scale factor 
    1616   USE eosbn2 , ONLY: eos                           ! eos routine 
    1717 
    18    USE lib_mpp, ONLY: ctl_stop ! ctl_stop routine 
    19    USE in_out_manager  !  
     18   USE lib_mpp, ONLY: ctl_stop                               ! ctl_stop routine 
     19   USE in_out_manager                                        !  
    2020 
    2121   IMPLICIT NONE 
     
    4545      ! ice shelf cavity 
    4646      SELECT CASE ( cn_isfload ) 
    47       CASE ( 'isomip' ) 
    48          CALL isf_load_isomip ( Kmm, pisfload ) 
     47      CASE ( 'uniform' ) 
     48         CALL isf_load_uniform ( Kmm, pisfload ) 
    4949      CASE DEFAULT 
    5050         CALL ctl_stop('STOP','method cn_isfload to compute ice shelf load does not exist (isomip), check your namelist') 
     
    5353   END SUBROUTINE isf_load 
    5454 
    55    SUBROUTINE isf_load_isomip( Kmm, pisfload ) 
     55   SUBROUTINE isf_load_uniform( Kmm, pisfload ) 
    5656      !!-------------------------------------------------------------------- 
    5757      !!                  ***  SUBROUTINE isf_load  *** 
     
    7979      znad = 1._wp                     !- To use density and not density anomaly 
    8080      ! 
    81       !                                !- assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 
    82       zts_top(:,:,jp_tem) = -1.9_wp   ;   zts_top(:,:,jp_sal) = 34.4_wp 
     81      !                                !- assume water displaced by the ice shelf is at T=rn_isfload_T and S=rn_isfload_S (rude) 
     82      zts_top(:,:,jp_tem) = rn_isfload_T   ;   zts_top(:,:,jp_sal) = rn_isfload_S 
    8383      ! 
    8484      DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf  
     
    113113      END DO 
    114114      ! 
    115    END SUBROUTINE isf_load_isomip 
     115   END SUBROUTINE isf_load_uniform 
    116116 
    117117END MODULE isfload 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfpar.F90

    r12068 r12077  
    1414   !!   isfpar       : compute ice shelf melt using a prametrisation of ice shelf cavities 
    1515   !!---------------------------------------------------------------------- 
    16    USE isf            ! ice shelf 
     16   USE isf_oce        ! ice shelf 
    1717   ! 
    1818   USE isfrst   , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine 
    1919   USE isftbl   , ONLY: isf_tbl_ktop, isf_tbl_lvl ! ice shelf top boundary layer properties subroutine 
     20   USE isfparmlt, ONLY: isfpar_mlt                ! ice shelf melt formulation subroutine 
     21   USE isfdiags , ONLY: isf_diags_flx             ! ice shelf diags subroutine 
    2022   USE isfutils , ONLY: debug, read_2dcstdta      ! ice shelf debug subroutine 
    21    USE isfparmlt, ONLY: isfpar_mlt     ! ice shelf melt formulation subroutine 
    22    USE isfdiags , ONLY: isf_diags_flx  ! ice shelf diags subroutine 
    2323   ! 
    2424   USE dom_oce  , ONLY: bathy          ! ocean space and time domain 
     
    4848      !! 
    4949      !! ** Purpose : compute the heat and fresh water due to ice shelf melting/freezing using a parametrisation  
     50      !! 
     51      !! ** Comment : in isf_par and all its call tree,  
     52      !!              'tbl' means parametrisation layer (ie how the far field temperature/salinity is computed)  
     53      !!              instead of in a proper top boundary layer as at the ice shelf ocean interface 
     54      !!              as the action to compute the properties of the tbl or the parametrisation layer are the same, 
     55      !!              (ie average T/S over a specific depth (can be across multiple levels)) 
     56      !!              the name tbl was kept. 
    5057      !! 
    5158      !!--------------------------------------------------------------------- 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfparmlt.F90

    r12068 r12077  
    88   !!---------------------------------------------------------------------- 
    99 
    10    USE isf            ! ice shelf 
     10   USE isf_oce                  ! ice shelf 
    1111   USE isftbl , ONLY: isf_tbl   ! ice shelf depth average 
    1212 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfstp.F90

    r12068 r12077  
    1414   !!---------------------------------------------------------------------- 
    1515   ! 
    16    USE isf            ! isf variables 
     16   USE isf_oce                                      ! isf variables 
    1717   USE isfload, ONLY: isf_load                      ! ice shelf load 
    1818   USE isftbl , ONLY: isf_tbl_lvl                   ! ice shelf boundary layer 
     
    187187            WRITE(numout,*) '      melt inside the cavity                  ln_isfcav_mlt   = ', ln_isfcav_mlt 
    188188            IF ( ln_isfcav_mlt) THEN 
    189                WRITE(numout,*) '         melt formulation                        cn_isfcav_mlt   = ', TRIM(cn_isfcav_mlt) 
    190                WRITE(numout,*) '         thickness of the top boundary layer     rn_htbl     = ', rn_htbl 
    191                WRITE(numout,*) '         gamma formulation                       cn_gammablk = ', TRIM(cn_gammablk)  
     189               WRITE(numout,*) '         melt formulation                         cn_isfcav_mlt= ', TRIM(cn_isfcav_mlt) 
     190               WRITE(numout,*) '         thickness of the top boundary layer      rn_htbl      = ', rn_htbl 
     191               WRITE(numout,*) '         gamma formulation                        cn_gammablk = ', TRIM(cn_gammablk)  
    192192               IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN  
    193                   WRITE(numout,*) '         gammat coefficient                       rn_gammat0  = ', rn_gammat0   
    194                   WRITE(numout,*) '         gammas coefficient                       rn_gammas0  = ', rn_gammas0   
    195                   WRITE(numout,*) '         top drag coef.    used (from namdrg_top) rn_Cd0      = ', r_Cdmin_top 
    196                   WRITE(numout,*) '         top background ke used (from namdrg_top) rn_ke0      = ', r_ke0_top 
     193                  WRITE(numout,*) '         gammat coefficient                       rn_gammat0   = ', rn_gammat0   
     194                  WRITE(numout,*) '         gammas coefficient                       rn_gammas0   = ', rn_gammas0   
     195                  WRITE(numout,*) '         top background ke used (from namdrg_top) rn_ke0       = ', r_ke0_top 
     196                  WRITE(numout,*) '         top drag coef.    used (from namdrg_top) rn_Cd0       = ', r_Cdmin_top 
    197197               END IF 
    198198            END IF 
     
    222222         END IF 
    223223 
    224          IF (ln_isfcav) WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload) 
     224         IF (ln_isfcav) THEN 
     225            WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload) 
     226            WRITE(numout,*) '         Temperature used to compute the ice shelf load            = ', rn_isfload_T 
     227            WRITE(numout,*) '         Salinity    used to compute the ice shelf load            = ', rn_isfload_S 
     228         END IF 
    225229         WRITE(numout,*) '' 
     230         FLUSH(numout) 
    226231 
    227232      END IF 
     
    245250      IF ( l_isfoasis .AND. ln_isf ) THEN 
    246251         ! 
    247          CALL ctl_stop( ' OASIS and ice shelf not tested' ) 
     252         CALL ctl_stop( ' ln_ctl and ice shelf not tested' ) 
    248253         ! 
    249254         ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation  
     
    277282      INTEGER               :: ios                  ! Local integer output status for namelist read 
    278283      !!---------------------------------------------------------------------- 
    279       NAMELIST/namisf/ ln_isf       ,                                                                               &  
    280          &             ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, sn_isfcav_fwf,   & 
    281          &             ln_isfpar_mlt, cn_isfpar_mlt, sn_isfpar_fwf, sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff, & 
    282          &             ln_isfcpl    , nn_drown, ln_isfcpl_cons, ln_isfdebug,                                        & 
    283          &             cn_isfload   , cn_isfdir 
     284      NAMELIST/namisf/ ln_isf        ,                                                           &  
     285         &             cn_gammablk   , rn_gammat0    , rn_gammas0    , rn_htbl, sn_isfcav_fwf,   & 
     286         &             ln_isfcav_mlt , cn_isfcav_mlt , sn_isfcav_fwf ,                           & 
     287         &             ln_isfpar_mlt , cn_isfpar_mlt , sn_isfpar_fwf ,                           & 
     288         &             sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff,                           & 
     289         &             ln_isfcpl     , nn_drown      , ln_isfcpl_cons, ln_isfdebug,              & 
     290         &             cn_isfload    , rn_isfload_T  , rn_isfload_S  , cn_isfdir 
    284291      !!---------------------------------------------------------------------- 
    285292      ! 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isftbl.F90

    r12068 r12077  
    1414   !!---------------------------------------------------------------------- 
    1515 
    16    USE isf    ! ice shelf variables 
     16   USE isf_oce ! ice shelf variables 
    1717 
    1818   USE dom_oce ! vertical scale factor and depth 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfutils.F90

    r11931 r12077  
    1515   USE lib_fortran   , ONLY: glob_sum, glob_min, glob_max             ! compute global value 
    1616   USE par_oce       , ONLY: jpi,jpj,jpk                              ! domain size 
    17    USE in_out_manager, ONLY: wp, lwp, numout                          ! miscelenious 
     17   USE dom_oce       , ONLY: nldi, nlei, nldj, nlej                   ! local domain 
     18   USE in_out_manager, ONLY: i8, wp, lwp, numout                          ! miscelenious 
     19   USE lib_mpp 
    1820 
    1921   IMPLICIT NONE 
     
    5456      !!                  ***  ROUTINE isf_debug2d  *** 
    5557      !! 
    56       !! ** Purpose : add debug print 
     58      !! ** Purpose : add debug print for 2d variables 
    5759      !! 
    5860      !!-------------------------- IN  ------------------------------------- 
    59       CHARACTER(LEN=256)          , INTENT(in   ) :: cdtxt 
     61      CHARACTER(LEN=*)            , INTENT(in   ) :: cdtxt 
    6062      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pvar 
    6163      !!-------------------------------------------------------------------- 
    62       REAL(wp) :: zmin, zmax, zsum 
     64      REAL(wp)    :: zmin, zmax, zsum 
     65      INTEGER(i8) :: imodd, ip 
     66      INTEGER     :: itmps,imods, ji, jj, jk 
    6367      !!-------------------------------------------------------------------- 
    6468      ! 
     69      ! global min/max/sum to check data range and NaN 
    6570      zsum = glob_sum( 'debug', pvar(:,:) ) 
    6671      zmin = glob_min( 'debug', pvar(:,:) ) 
    6772      zmax = glob_max( 'debug', pvar(:,:) ) 
    6873      ! 
    69       IF (lwp) WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum 
     74      ! basic check sum to check reproducibility 
     75      ! TRANSFER function find out the integer corresponding to pvar(i,j) bit pattern 
     76      ! MOD allow us to keep only the latest digits during the sum 
     77      ! imod is not choosen to be very large as at the end there is a classic mpp_sum 
     78      imodd=65521 ! highest prime number < 2**16 with i8 type 
     79      imods=65521 ! highest prime number < 2**16 with default integer for mpp_sum subroutine 
     80      itmps=0 
     81      DO jj=nldj,nlej 
     82         DO ji=nldi,nlei 
     83            itmps = MOD(itmps + MOD(TRANSFER(pvar(ji,jj), ip),imodd), imods) 
     84         END DO 
     85      END DO 
     86      CALL mpp_sum('debug',itmps) 
    7087      ! 
    71       CALL FLUSH(numout) 
     88      ! print out 
     89      IF (lwp) THEN 
     90         WRITE(numout,*) TRIM(cdtxt),' (min, max, sum, tag) : ',zmin, zmax, zsum, itmps 
     91         CALL FLUSH(numout) 
     92      END IF 
    7293      ! 
    7394   END SUBROUTINE debug2d 
     
    7798      !!                  ***  ROUTINE isf_debug3d  *** 
    7899      !! 
    79       !! ** Purpose : add debug print 
     100      !! ** Purpose : add debug print for 3d variables 
    80101      !! 
    81102      !!-------------------------- IN  ------------------------------------- 
    82       CHARACTER(LEN=256)              , INTENT(in   ) :: cdtxt 
     103      CHARACTER(LEN=*)                , INTENT(in   ) :: cdtxt 
    83104      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pvar 
    84105      !!-------------------------------------------------------------------- 
    85106      REAL(wp) :: zmin, zmax, zsum 
     107      INTEGER(i8) :: imodd, ip 
     108      INTEGER     :: itmps,imods, ji, jj, jk 
    86109      !!-------------------------------------------------------------------- 
    87110      ! 
     111      ! global min/max/sum to check data range and NaN 
    88112      zsum = glob_sum( 'debug', pvar(:,:,:) ) 
    89113      zmin = glob_min( 'debug', pvar(:,:,:) ) 
    90114      zmax = glob_max( 'debug', pvar(:,:,:) ) 
    91115      ! 
    92       IF (lwp) WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum 
     116      ! basic check sum to check reproducibility 
     117      ! TRANSFER function find out the integer corresponding to pvar(i,j) bit pattern 
     118      ! MOD allow us to keep only the latest digits during the sum 
     119      ! imod is not choosen to be very large as at the end there is a classic mpp_sum 
     120      imodd=65521 ! highest prime number < 2**16 with i8 type 
     121      imods=65521 ! highest prime number < 2**16 with default integer for mpp_sum subroutine 
     122      itmps=0 
     123      DO jk=1,jpk 
     124         DO jj=nldj,nlej 
     125            DO ji=nldi,nlei 
     126               itmps = MOD(itmps + MOD(TRANSFER(pvar(ji,jj,jk), ip),imodd), imods) 
     127            END DO 
     128         END DO 
     129      END DO 
     130      CALL mpp_sum('debug',itmps) 
    93131      ! 
    94       CALL FLUSH(numout) 
     132      ! print out 
     133      IF (lwp) THEN 
     134         WRITE(numout,*) TRIM(cdtxt),' (min, max, sum, tag) : ',zmin, zmax, zsum, itmps 
     135         CALL FLUSH(numout) 
     136      END IF 
    95137      ! 
    96138   END SUBROUTINE debug3d 
    97139 
    98140END MODULE isfutils 
    99  
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldfslp.F90

    r12068 r12077  
    2121   !!---------------------------------------------------------------------- 
    2222   USE oce            ! ocean dynamics and tracers 
    23    USE isf            ! ice shelf 
     23   USE isf_oce        ! ice shelf 
    2424   USE dom_oce        ! ocean space and time domain 
    2525!   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldftra.F90

    r11822 r12077  
    664664                  ! eddies using the isopycnal slopes calculated in ldfslp.F :  
    665665                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    666                   ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
     666                  ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
    667667                  zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 
    668668                  zhw(ji,jj) = zhw(ji,jj) + ze3w 
     
    682682                  ! eddies using the isopycnal slopes calculated in ldfslp.F :  
    683683                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    684                   ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
     684                  ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
    685685                  zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    686686                     &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbccpl.F90

    r12068 r12077  
    2727   USE sbcwave         ! surface boundary condition: waves 
    2828   USE phycst          ! physical constants 
     29   USE isf_oce , ONLY : l_isfoasis, fwfisf_oasis ! ice shelf boundary condition 
    2930#if defined key_si3 
    3031   USE ice            ! ice variables 
     
    3637   USE eosbn2         !  
    3738   USE sbcrnf  , ONLY : l_rnfcpl 
    38    USE isf     , ONLY : ln_isf, l_isfoasis, fwfisf_oasis 
    3939#if defined key_cice 
    4040   USE ice_domain_size, only: ncat 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcfwb.F90

    r12068 r12077  
    1717   USE dom_oce        ! ocean space and time domain 
    1818   USE sbc_oce        ! surface ocean boundary condition 
     19   USE isf_oce , ONLY : fwfisf_cav, fwfisf_par                    ! ice shelf melting contribution 
    1920   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass 
    2021   USE phycst         ! physical constants 
    2122   USE sbcrnf         ! ocean runoffs 
    22    USE isf            ! ice shelf melting contribution 
    2323   USE sbcssr         ! Sea-Surface damping terms 
    2424   ! 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcrnf.F90

    r12068 r12077  
    1919   USE phycst         ! physical constants 
    2020   USE sbc_oce        ! surface boundary condition variables 
    21    USE isf            ! ice shelf 
    2221   USE eosbn2         ! Equation Of State 
    2322   USE closea         ! closed seas 
     
    127126               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
    128127            END WHERE 
    129             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg 
    130                rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rLfusisf * r1_rau0_rcp 
    131             END WHERE 
    132128         ELSE                                                        ! use SST as runoffs temperature 
    133129            !CEOD River is fresh water so must at least be 0 unless we consider ice 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcssm.F90

    r11027 r12077  
    6161      ! 
    6262      !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity) 
    63       DO jj = 1, jpj 
    64          DO ji = 1, jpi 
    65             zts(ji,jj,jp_tem) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) 
    66             zts(ji,jj,jp_sal) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) 
    67          END DO 
    68       END DO 
     63      zts(:,:,jp_tem) = ts(:,:,1,jp_tem,Kmm) 
     64      zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm) 
    6965      ! 
    7066      IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        ! 
     
    7369         ssv_m(:,:) = vv(:,:,1,Kbb) 
    7470         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) ) 
    75          ELSE                    ;   sst_m(:,:) = zts(:,:,jp_tem) 
     71         ELSE                   ;   sst_m(:,:) = zts(:,:,jp_tem) 
    7672         ENDIF 
    7773         sss_m(:,:) = zts(:,:,jp_sal) 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/eosbn2.F90

    r11822 r12077  
    2929   !!   eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass 
    3030   !!   eos_insitu_2d : Compute the in situ density for 2d fields 
    31    !!   bn2           : Compute the Brunt-Vaisala frequency 
    3231   !!   bn2           : compute the Brunt-Vaisala frequency 
    3332   !!   eos_pt_from_ct: compute the potential temperature from the Conservative Temperature 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traatf.F90

    r12068 r12077  
    2929   USE sbc_oce         ! surface boundary condition: ocean 
    3030   USE sbcrnf          ! river runoffs 
    31    USE isf             ! ice shelf melting 
     31   USE isf_oce         ! ice shelf melting 
    3232   USE zdf_oce         ! ocean vertical mixing 
    3333   USE domvvl          ! variable volume 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traisf.F90

    r12068 r12077  
    1010   !!   tra_isf       : update the tracer trend at ocean surface 
    1111   !!---------------------------------------------------------------------- 
    12    USE oce            ! ocean dynamics and active tracers 
    13    USE dom_oce        ! ocean space domain variables 
    14    USE phycst         ! physical constant 
    15    USE eosbn2         ! Equation Of State 
    16    USE isf            ! Ice shelf variable 
    17    USE isfutils       ! 
    18    ! 
    19    USE in_out_manager ! I/O manager 
    20    USE iom            ! xIOS server 
    21    USE timing         ! Timing 
     12   USE isf_oce                                     ! Ice shelf variables 
     13   USE dom_oce , ONLY : e3t, r1_e1e2t            ! ocean space domain variables 
     14   USE isfutils, ONLY : debug                      ! debug option 
     15   USE timing  , ONLY : timing_start, timing_stop  ! Timing 
     16   USE in_out_manager                              ! I/O manager 
    2217 
    2318   IMPLICIT NONE 
     
    6560      IF ( ln_isfcpl ) THEN 
    6661         ! 
     62         ! Dynamical stability at start up after change in under ice shelf cavity geometry is achieve by correcting the divergence. 
     63         ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping  
     64         ! the same as at the end of the latest time step. So correction need to be apply at nit000 (euler time step) and 
     65         ! half of it at nit000+1 (leap frog time step). 
     66         ! in accordance to this, the heat content flux due to injected water need to be added in the temperature and salt trend 
     67         ! at time step nit000 and nit000+1 
    6768         IF ( kt == nit000  ) CALL tra_isf_cpl(Kmm, risfcpl_tsc       , pts(:,:,:,:,Krhs)) 
    6869         IF ( kt == nit000+1) CALL tra_isf_cpl(Kmm, risfcpl_tsc*0.5_wp, pts(:,:,:,:,Krhs)) 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ZDF/zdfmxl.F90

    r12068 r12077  
    1212   !!---------------------------------------------------------------------- 
    1313   USE oce            ! ocean dynamics and tracers variables 
    14    USE isf            ! ice shelf 
     14   USE isf_oce        ! ice shelf 
    1515   USE dom_oce        ! ocean space and time domain variables 
    1616   USE trc_oce  , ONLY: l_offline         ! ocean space and time domain variables 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step.F90

    r12068 r12077  
    115115 
    116116      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    117       ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice) 
    118       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    119       IF( ln_isf     )   CALL isf_stp ( kstp, Nnn ) 
     117      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice) 
     118      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    120119      IF( ln_tide    )   CALL sbc_tide( kstp )                        ! update tide potential 
    121120      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)  
    122121      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn, kt_offset = +1 )   ! update dynamic & tracer data at open boundaries 
     122      IF( ln_isf     )   CALL isf_stp ( kstp, Nnn ) 
    123123                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition (including sea-ice) 
    124124 
  • NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step_oce.F90

    r12068 r12077  
    2222   USE sbcwave         ! Wave intialisation 
    2323 
    24    USE isf 
     24   USE isf_oce         ! ice shelf boundary condition 
    2525   USE isfstp          ! ice shelf boundary condition     (isf_stp routine) 
    2626 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/EXPREF/namelist_cfg

    r11896 r12077  
    165165   ! ---------------- ice shelf load ------------------------------- 
    166166   ! 
    167    cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc) 
     167   cn_isfload = 'uniform'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc) 
     168      rn_isfload_T = -1.0 
     169      rn_isfload_S =  34.2 
    168170   ! 
    169171   ! ---------------- ice shelf melt formulation ------------------------------- 
     
    182184         !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf 
    183185         !              !  cn_isfcav_mlt = 2eq or 3eq cases: 
    184          cn_gammablk = 'ad15'    ! scheme to compute gammat/s (spe,ad15,hj99) 
    185          !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
    186          !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999) 
     186         cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99) 
     187         !                       ! spe      = constant transfert velocity (rn_gammat0, rn_gammas0) 
     188         !                       ! vel      = velocity dependent transfert velocity (u* * gammat/s) (Asay-Davis et al. 2016 for a short description) 
     189         !                       ! vel_stab = velocity and stability dependent transfert coeficient (Holland et al. 1999 for a complete description) 
    187190         rn_gammat0  = 0.0215    ! gammat coefficient used in blk formula 
    188191         rn_gammas0  = 0.614e-3  ! gammas coefficient used in blk formula 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/dtatsd.F90

    r11889 r12077  
    6969      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
    7070      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
    71 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist', lwp ) 
     71901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist' ) 
    7272      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
    7373      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
    74 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp ) 
     74902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist' ) 
    7575      IF(lwm) WRITE ( numond, namtsd ) 
    7676 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/eosbn2.F90

    r11889 r12077  
    2929   !!   eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass 
    3030   !!   eos_insitu_2d : Compute the in situ density for 2d fields 
    31    !!   bn2           : Compute the Brunt-Vaisala frequency 
     31   !!   bn2           : compute the Brunt-Vaisala frequency 
     32   !!   eos_pt_from_ct: compute the potential temperature from the Conservative Temperature 
    3233   !!   eos_rab       : generic interface of in situ thermal/haline expansion ratio  
    3334   !!   eos_rab_3d    : compute in situ thermal/haline expansion ratio 
     
    7475 
    7576   !                               !!** Namelist nameos ** 
    76    LOGICAL , PUBLIC ::   ln_TEOS10   ! determine if eos_pt_from_ct is used to compute sst_m 
    77    LOGICAL , PUBLIC ::   ln_EOS80   ! determine if eos_pt_from_ct is used to compute sst_m 
    78    LOGICAL , PUBLIC ::   ln_SEOS   ! determine if eos_pt_from_ct is used to compute sst_m 
     77   LOGICAL , PUBLIC ::   ln_TEOS10 
     78   LOGICAL , PUBLIC ::   ln_EOS80 
     79   LOGICAL , PUBLIC ::   ln_SEOS 
    7980   LOGICAL , PUBLIC ::   ln_LEOS   ! determine if linear eos is used  
    8081 
     
    624625 
    625626 
    626    SUBROUTINE rab_3d( pts, pab ) 
     627   SUBROUTINE rab_3d( pts, pab, Kmm ) 
    627628      !!---------------------------------------------------------------------- 
    628629      !!                 ***  ROUTINE rab_3d  *** 
     
    634635      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
    635636      !!---------------------------------------------------------------------- 
     637      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    636638      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity 
    637639      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio 
     
    652654               DO ji = 1, jpi 
    653655                  ! 
    654                   zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth 
     656                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
    655657                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    656658                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     
    710712                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    711713                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    712                   zh  = gdept_n(ji,jj,jk)                ! depth in meters at t-point 
     714                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point 
    713715                  ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
    714716                  ! 
     
    730732                  zt  = pts (ji,jj,jk,jp_tem) - (-1._wp) 
    731733                  zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    732                   zh  = gdept_n(ji,jj,jk)                 ! depth in meters at t-point 
     734                  zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point 
    733735                  ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask 
    734736                  ! 
     
    757759 
    758760 
    759    SUBROUTINE rab_2d( pts, pdep, pab ) 
     761   SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 
    760762      !!---------------------------------------------------------------------- 
    761763      !!                 ***  ROUTINE rab_2d  *** 
     
    765767      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
    766768      !!---------------------------------------------------------------------- 
     769      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    767770      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity 
    768771      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m] 
     
    891894 
    892895 
    893    SUBROUTINE rab_0d( pts, pdep, pab ) 
     896   SUBROUTINE rab_0d( pts, pdep, pab, Kmm ) 
    894897      !!---------------------------------------------------------------------- 
    895898      !!                 ***  ROUTINE rab_0d  *** 
     
    899902      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
    900903      !!---------------------------------------------------------------------- 
     904      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    901905      REAL(wp), DIMENSION(jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity 
    902906      REAL(wp),                      INTENT(in   ) ::   pdep   ! depth                  [m] 
     
    9991003 
    10001004 
    1001    SUBROUTINE bn2( pts, pab, pn2 ) 
     1005   SUBROUTINE bn2( pts, pab, pn2, Kmm ) 
    10021006      !!---------------------------------------------------------------------- 
    10031007      !!                  ***  ROUTINE bn2  *** 
     
    10131017      !! 
    10141018      !!---------------------------------------------------------------------- 
     1019      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    10151020      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu] 
    10161021      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1] 
     
    10261031         DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90 
    10271032            DO ji = 1, jpi 
    1028                zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
    1029                   &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )  
     1033               zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   & 
     1034                  &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )  
    10301035                  ! 
    10311036               zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
     
    10341039               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    10351040                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    1036                   &            / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
     1041                  &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
    10371042            END DO 
    10381043         END DO 
     
    12031208 
    12041209 
    1205    SUBROUTINE eos_pen( pts, pab_pe, ppen ) 
     1210   SUBROUTINE eos_pen( pts, pab_pe, ppen, Kmm ) 
    12061211      !!---------------------------------------------------------------------- 
    12071212      !!                 ***  ROUTINE eos_pen  *** 
     
    12231228      !!                    pab_pe(:,:,:,jp_sal) is beta_pe 
    12241229      !!---------------------------------------------------------------------- 
     1230      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    12251231      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity 
    12261232      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe 
     
    12421248               DO ji = 1, jpi 
    12431249                  ! 
    1244                   zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth 
     1250                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
    12451251                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    12461252                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     
    13061312                  zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
    13071313                  zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
    1308                   zh  = gdept_n(ji,jj,jk)              ! depth in meters  at t-point 
     1314                  zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point 
    13091315                  ztm = tmask(ji,jj,jk)                ! tmask 
    13101316                  zn  = 0.5_wp * zh * r1_rau0 * ztm 
     
    13261332                  zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0) 
    13271333                  zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    1328                   zh  = gdept_n(ji,jj,jk)                ! depth in meters  at t-point 
     1334                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point 
    13291335                  ztm = tmask(ji,jj,jk)                  ! tmask 
    13301336                  zn  = 0.5_wp * zh * r1_rau0 * ztm 
     
    13681374      REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state 
    13691375      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 
    1370 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist', lwp ) 
     1376901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' ) 
    13711377      ! 
    13721378      REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state 
    13731379      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 
    1374 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp ) 
     1380902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' ) 
    13751381      IF(lwm) WRITE( numond, nameos ) 
    13761382      ! 
     
    17821788         ! 
    17831789      CASE( np_seos )                        !==  Simplified EOS     ==! 
     1790 
     1791         r1_S0  = 0.875_wp/35.16504_wp   ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct) 
     1792          
    17841793         IF(lwp) THEN 
    17851794            WRITE(numout,*) 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfcavgam.F90

    r11889 r12077  
    1010   !!   isfcav_gammats       : compute exchange coeficient gamma  
    1111   !!---------------------------------------------------------------------- 
    12    USE isf 
     12   USE isf_oce 
    1313   USE isfutils, ONLY: debug 
    1414   USE isftbl  , ONLY: isf_tbl 
    1515 
    16    USE oce     , ONLY: un, vn, rn2         ! ocean dynamics and tracers 
     16   USE oce     , ONLY: uu, vv, rn2         ! ocean dynamics and tracers 
    1717   USE phycst  , ONLY: grav, vkarmn        ! physical constant 
    1818   USE eosbn2  , ONLY: eos_rab             ! equation of state 
     
    4141   !!----------------------------------------------------------------------------------------------------- 
    4242   ! 
    43    SUBROUTINE isfcav_gammats( pttbl, pstbl, pqoce, pqfwf, pgt, pgs ) 
     43   SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pgt, pgs ) 
    4444      !!---------------------------------------------------------------------- 
    4545      !! ** Purpose    : compute the coefficient echange for heat and fwf flux 
    4646      !! 
    4747      !! ** Method     : select the gamma formulation 
    48       !!                 3 method available (cst, AD15 and HJ99) 
     48      !!                 3 method available (cst, vel and vel_stab) 
    4949      !!--------------------------------------------------------------------- 
    5050      !!-------------------------- OUT ------------------------------------- 
    5151      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt  , pgs      ! gamma t and gamma s  
    5252      !!-------------------------- IN  ------------------------------------- 
     53      INTEGER                                     :: Kmm             ! ocean time level index 
    5354      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf    ! isf heat and fwf 
    5455      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl    ! top boundary layer tracer 
     
    5758      !!--------------------------------------------------------------------- 
    5859      ! 
    59       ! compute velocity in the tbl if needed 
     60      !========================================== 
     61      ! 1.: compute velocity in the tbl if needed 
     62      !========================================== 
     63      ! 
    6064      SELECT CASE ( cn_gammablk ) 
    6165      CASE ( 'spe'  )  
    6266         ! gamma is constant (specified in namelist) 
    6367         ! nothing to do 
    64       CASE ('ad15', 'hj99') 
     68      CASE ('vel', 'vel_stab') 
    6569         ! compute velocity in tbl 
    66          CALL isf_tbl(un(:,:,:) ,zutbl(:,:),'U', miku, rhisf_tbl_cav) 
    67          CALL isf_tbl(vn(:,:,:) ,zvtbl(:,:),'V', mikv, rhisf_tbl_cav) 
     70         CALL isf_tbl(Kmm, uu(:,:,:,Kmm) ,zutbl(:,:),'U', miku, rhisf_tbl_cav) 
     71         CALL isf_tbl(Kmm, vv(:,:,:,Kmm) ,zvtbl(:,:),'V', mikv, rhisf_tbl_cav) 
    6872         ! 
    6973         ! mask velocity in tbl with ice shelf mask 
     
    7882      END SELECT 
    7983      !  
    80       ! compute gamma 
     84      !========================================== 
     85      ! 2.: compute gamma 
     86      !========================================== 
     87      ! 
    8188      SELECT CASE ( cn_gammablk ) 
    8289      CASE ( 'spe'  ) ! gamma is constant (specified in namelist) 
    8390         pgt(:,:) = rn_gammat0 
    8491         pgs(:,:) = rn_gammas0 
    85       CASE ( 'ad15' ) ! gamma is proportional to u* 
    86          CALL gammats_AD15 (              zutbl, zvtbl, rCd0_top, rn_vtide**2,               pgt, pgs ) 
    87       CASE ( 'hj99' ) ! gamma depends of stability of boundary layer and u* 
    88          CALL gammats_HJ99 (pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top  , pqoce, pqfwf, pgt, pgs ) 
     92      CASE ( 'vel' ) ! gamma is proportional to u* 
     93         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs ) 
     94      CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u* 
     95         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs ) 
    8996      CASE DEFAULT 
    9097         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)') 
    9198      END SELECT 
    9299      ! 
    93       ! ouput exchange coeficient and tbl velocity 
     100      !========================================== 
     101      ! 3.: output and debug 
     102      !========================================== 
     103      ! 
    94104      CALL iom_put('isfgammat', pgt(:,:)) 
    95105      CALL iom_put('isfgammas', pgs(:,:)) 
     
    106116   !!----------------------------------------------------------------------------------------------------- 
    107117   ! 
    108    SUBROUTINE gammats_AD15(putbl, pvtbl, pCd, pke2, &   ! <<== in 
     118   SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, &   ! <<== in 
    109119      &                                  pgt, pgs   )   ! ==>> out gammats [m/s] 
    110120      !!---------------------------------------------------------------------- 
     
    113123      !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar ) 
    114124      !! 
    115       !! ** Reference  : Jenkins et al., 2010, JPO, p2298-2312 
    116       !!                 Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016 
     125      !! ** Reference  : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016 
    117126      !!--------------------------------------------------------------------- 
    118127      !!-------------------------- OUT ------------------------------------- 
     
    136145      CALL iom_put('isfustar',zustar(:,:)) 
    137146      ! 
    138    END SUBROUTINE gammats_AD15 
    139  
    140    SUBROUTINE gammats_HJ99( pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in 
    141       &                                                            pgt  , pgs    )  ! ==>> out gammats [m/s] 
     147   END SUBROUTINE gammats_vel 
     148 
     149   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in 
     150      &                                                                     pgt  , pgs    )  ! ==>> out gammats [m/s] 
    142151      !!---------------------------------------------------------------------- 
    143152      !! ** Purpose    : compute the coefficient echange coefficient  
     
    150159      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas 
    151160      !!-------------------------- IN  ------------------------------------- 
     161      INTEGER                                     :: Kmm            ! ocean time level index 
    152162      REAL(wp),                     INTENT(in   ) :: pke2           ! background velocity squared 
    153163      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf   ! surface heat flux and fwf flux 
     
    199209               ! compute Rc number (as done in zdfric.F90) 
    200210!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation 
    201                zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1) 
     211               zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm) 
    202212               !                                            ! shear of horizontal velocity 
    203                zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  & 
    204                   &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
    205                zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  & 
    206                   &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
     213               zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  & 
     214                  &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  ) 
     215               zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  & 
     216                  &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  ) 
    207217               !                                            ! richardson number (minimum value set to zero) 
    208218               zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 
     
    211221               zts(jp_tem) = pttbl(ji,jj) 
    212222               zts(jp_sal) = pstbl(ji,jj) 
    213                zdep        = gdepw_n(ji,jj,ikt) 
    214                ! 
    215                CALL eos_rab( zts, zdep, zab ) 
     223               zdep        = gdepw(ji,jj,ikt,Kmm) 
     224               ! 
     225               CALL eos_rab( zts, zdep, zab, Kmm ) 
    216226               ! 
    217227               ! compute length scale (Eq ??) 
     
    220230               ! compute Monin Obukov Length 
    221231               ! Maximum boundary layer depth (Eq ??) 
    222                zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp 
     232               zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp 
    223233               ! 
    224234               ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??) 
     
    243253      END DO 
    244254 
    245    END SUBROUTINE gammats_HJ99 
     255   END SUBROUTINE gammats_vel_stab 
    246256 
    247257END MODULE isfcavgam 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfstp.F90

    r11931 r12077  
    1414   !!---------------------------------------------------------------------- 
    1515   ! 
    16    USE isf            ! isf variables 
     16   USE isf_oce                                      ! isf variables 
    1717   USE isfload, ONLY: isf_load                      ! ice shelf load 
    1818   USE isftbl , ONLY: isf_tbl_lvl                   ! ice shelf boundary layer 
     
    2121   USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables 
    2222 
    23    USE dom_oce, ONLY: ht_n, e3t_n, ln_isfcav, ln_linssh ! ocean space and time domain 
     23   USE dom_oce, ONLY: ht, e3t, ln_isfcav, ln_linssh    ! ocean space and time domain 
    2424   USE domvvl,  ONLY: ln_vvl_zstar                      ! zstar logical 
    2525   USE zdfdrg,  ONLY: r_Cdmin_top, r_ke0_top            ! vertical physics: top/bottom drag coef. 
    2626   ! 
    2727   USE lib_mpp, ONLY: ctl_stop, ctl_nam 
     28   USE fldread, ONLY: FLD, FLD_N 
    2829   USE in_out_manager ! I/O manager 
    2930   USE timing 
     
    4243CONTAINS 
    4344  
    44   SUBROUTINE isf_stp( kt ) 
     45  SUBROUTINE isf_stp( kt, Kmm ) 
    4546      !!--------------------------------------------------------------------- 
    4647      !!                  ***  ROUTINE isf_stp  *** 
     
    5859      !!---------------------------------------------------------------------- 
    5960      INTEGER, INTENT(in) ::   kt   ! ocean time step 
     61      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    6062      !!--------------------------------------------------------------------- 
    6163      ! 
    6264      IF( ln_timing )   CALL timing_start('isf') 
    6365      ! 
     66      !======================================================================= 
     67      ! 1.: compute melt and associated heat fluxes in the ice shelf cavities 
     68      !======================================================================= 
     69      ! 
    6470      IF ( ln_isfcav_mlt ) THEN 
    6571         ! 
    66          ! before time step  
     72         ! 1.1: before time step  
    6773         IF ( kt /= nit000 ) THEN  
    6874            risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:) 
     
    7076         END IF 
    7177         ! 
    72          ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 
     78         ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 
    7379         rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:) 
    74          CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 
    75          ! 
    76          ! compute ice shelf melt 
    77          CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav) 
     80         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 
     81         ! 
     82         ! 1.3: compute ice shelf melt 
     83         CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav) 
    7884         ! 
    7985      END IF 
    8086      !  
     87      !================================================================================= 
     88      ! 2.: compute melt and associated heat fluxes for not resolved ice shelf cavities 
     89      !================================================================================= 
     90      ! 
    8191      IF ( ln_isfpar_mlt ) THEN 
    8292         ! 
    83          ! before time step  
     93         ! 2.1: before time step  
    8494         IF ( kt /= nit000 ) THEN  
    8595            risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:) 
     
    8797         END IF 
    8898         ! 
    89          ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 
    90          ! by simplicity, we assume the top level where param applied do not change with time 
     99         ! 2.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 
     100         ! by simplicity, we assume the top level where param applied do not change with time (done in init part) 
    91101         rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:) 
    92          CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 
    93          ! 
    94          ! compute ice shelf melt 
    95          CALL isf_par( kt, risf_par_tsc, fwfisf_par) 
    96          ! 
    97       END IF 
    98       ! 
    99       IF ( ln_isfcpl ) THEN 
    100          ! after step nit000 + 2 we do not need anymore the risfcpl_ arrays 
    101          IF ( kt == nit000 + 2 ) CALL isf_dealloc_cpl() 
    102  
    103          IF ( lrst_oce ) CALL isfcpl_rst_write(kt) 
    104       END IF 
     102         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 
     103         ! 
     104         ! 2.3: compute ice shelf melt 
     105         CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par) 
     106         ! 
     107      END IF 
     108      ! 
     109      !================================================================================== 
     110      ! 3.: output specific restart variable in case of coupling with an ice sheet model 
     111      !================================================================================== 
     112      ! 
     113      IF ( ln_isfcpl .AND. lrst_oce ) CALL isfcpl_rst_write(kt, Kmm) 
    105114      ! 
    106115      IF( ln_timing )   CALL timing_stop('isf') 
     
    108117   END SUBROUTINE isf_stp 
    109118 
    110    SUBROUTINE isf_init 
     119   SUBROUTINE isf_init(Kbb, Kmm, Kaa) 
    111120      !!--------------------------------------------------------------------- 
    112121      !!                  ***  ROUTINE isfstp_init  *** 
     
    122131      !!              - call cav/param/isfcpl init routine 
    123132      !!---------------------------------------------------------------------- 
     133      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa      ! ocean time level indices 
    124134      ! 
    125135      ! constrain: l_isfoasis need to be known 
     
    135145      ! 
    136146      ! compute ice shelf load 
    137       IF ( ln_isfcav ) CALL isf_load( risfload ) 
     147      IF ( ln_isfcav ) CALL isf_load( Kmm, risfload ) 
    138148      ! 
    139149      ! terminate routine now if no ice shelf melt formulation specify 
     
    150160         !--------------------------------------------------------------------------------------------------------------------- 
    151161         ! initialisation ice sheet coupling 
    152          IF ( ln_isfcpl ) CALL isfcpl_init() 
     162         IF( ln_isfcpl ) CALL isfcpl_init(Kbb, Kmm, Kaa) 
    153163         ! 
    154164      END IF 
     
    177187            WRITE(numout,*) '      melt inside the cavity                  ln_isfcav_mlt   = ', ln_isfcav_mlt 
    178188            IF ( ln_isfcav_mlt) THEN 
    179                WRITE(numout,*) '         melt formulation                        cn_isfcav_mlt   = ', TRIM(cn_isfcav_mlt) 
    180                WRITE(numout,*) '         thickness of the top boundary layer     rn_htbl     = ', rn_htbl 
    181                WRITE(numout,*) '         gamma formulation                       cn_gammablk = ', TRIM(cn_gammablk)  
     189               WRITE(numout,*) '         melt formulation                         cn_isfcav_mlt= ', TRIM(cn_isfcav_mlt) 
     190               WRITE(numout,*) '         thickness of the top boundary layer      rn_htbl      = ', rn_htbl 
     191               WRITE(numout,*) '         gamma formulation                        cn_gammablk = ', TRIM(cn_gammablk)  
    182192               IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN  
    183                   WRITE(numout,*) '         gammat coefficient                       rn_gammat0  = ', rn_gammat0   
    184                   WRITE(numout,*) '         gammas coefficient                       rn_gammas0  = ', rn_gammas0   
    185                   WRITE(numout,*) '         top background ke used (from namdrg_top) rn_vtide**2 = ', rn_vtide**2 
    186                   WRITE(numout,*) '         top drag coef.    used (from namdrg_top) rn_Cd0      = ', r_Cdmin_top 
     193                  WRITE(numout,*) '         gammat coefficient                       rn_gammat0   = ', rn_gammat0   
     194                  WRITE(numout,*) '         gammas coefficient                       rn_gammas0   = ', rn_gammas0   
     195                  WRITE(numout,*) '         top background ke used (from namdrg_top) rn_vtide**2  = ', rn_vtide**2 
     196                  WRITE(numout,*) '         top drag coef.    used (from namdrg_top) rn_Cd0       = ', r_Cdmin_top 
    187197               END IF 
    188198            END IF 
     
    212222         END IF 
    213223 
    214          IF (ln_isfcav) WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload) 
     224         IF (ln_isfcav) THEN 
     225            WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload) 
     226            WRITE(numout,*) '         Temperature used to compute the ice shelf load            = ', rn_isfload_T 
     227            WRITE(numout,*) '         Salinity    used to compute the ice shelf load            = ', rn_isfload_S 
     228         END IF 
    215229         WRITE(numout,*) '' 
     230         FLUSH(numout) 
    216231 
    217232      END IF 
     
    267282      INTEGER               :: ios                  ! Local integer output status for namelist read 
    268283      !!---------------------------------------------------------------------- 
    269       NAMELIST/namisf/ ln_isf       ,                                                                               &  
    270          &             ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, sn_isfcav_fwf,   & 
    271          &             ln_isfpar_mlt, cn_isfpar_mlt, sn_isfpar_fwf, sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff, & 
    272          &             ln_isfcpl    , nn_drown, ln_isfcpl_cons, ln_isfdebug, rn_vtide,                              & 
    273          &             cn_isfload   , cn_isfdir 
     284      NAMELIST/namisf/ ln_isf        ,                                                           &  
     285         &             cn_gammablk   , rn_gammat0    , rn_gammas0    , rn_htbl, sn_isfcav_fwf,   & 
     286         &             ln_isfcav_mlt , cn_isfcav_mlt , sn_isfcav_fwf ,                           & 
     287         &             ln_isfpar_mlt , cn_isfpar_mlt , sn_isfpar_fwf ,                           & 
     288         &             sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff,                           & 
     289         &             ln_isfcpl     , nn_drown      , ln_isfcpl_cons, ln_isfdebug, rn_vtide,    & 
     290         &             cn_isfload    , rn_isfload_T  , rn_isfload_S  , cn_isfdir 
    274291      !!---------------------------------------------------------------------- 
    275292      ! 
    276293      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
    277294      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901) 
    278 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist', lwp ) 
     295901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist' ) 
    279296      ! 
    280297      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
    281298      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 ) 
    282 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp ) 
     299902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist' ) 
    283300      IF(lwm) WRITE ( numond, namisf ) 
    284301 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/istate.F90

    r11889 r12077  
    5050CONTAINS 
    5151 
    52    SUBROUTINE istate_init 
     52   SUBROUTINE istate_init( Kbb, Kmm, Kaa ) 
    5353      !!---------------------------------------------------------------------- 
    5454      !!                   ***  ROUTINE istate_init  *** 
     
    5656      !! ** Purpose :   Initialization of the dynamics and tracer fields. 
    5757      !!---------------------------------------------------------------------- 
     58      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices 
     59      ! 
    5860      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    5961!!gm see comment further down 
     
    7577      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    7678      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
    77       tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk 
     79      ts   (:,:,:,:,Kaa) = 0._wp                               ! set one for all to 0 at level jpk 
    7880      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    7981#if defined key_agrif 
    80       ua   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization 
    81       va   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization     
     82      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization 
     83      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization     
    8284#endif 
    8385 
    8486      IF( ln_rstart ) THEN                    ! Restart from a file 
    8587         !                                    ! ------------------- 
    86          CALL rst_read                        ! Read the restart file 
     88         CALL rst_read( Kbb, Kmm )            ! Read the restart file 
    8789         CALL day_init                        ! model calendar (using both namelist and restart infos) 
    8890         ! 
     
    9597         ! 
    9698         IF( ln_tsd_init ) THEN                
    97             CALL dta_tsd( nit000, 'ini', tsb )       ! read 3D T and S data at nit000 
     99            CALL dta_tsd( nit000, 'ini', ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000 
    98100            ! 
    99             sshb(:,:)   = 0._wp               ! set the ocean at rest 
     101            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest 
    100102            IF( ll_wd ) THEN 
    101                sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
     103               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
    102104               ! 
    103105               ! Apply minimum wetdepth criterion 
     
    105107               DO jj = 1,jpj 
    106108                  DO ji = 1,jpi 
    107                      IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN 
    108                         sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
     109                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
     110                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
    109111                     ENDIF 
    110112                  END DO 
    111113               END DO  
    112114            ENDIF  
    113             ub  (:,:,:) = 0._wp 
    114             vb  (:,:,:) = 0._wp   
     115            uu  (:,:,:,Kbb) = 0._wp 
     116            vv  (:,:,:,Kbb) = 0._wp   
    115117            ! 
    116118         ELSE                                 ! user defined initial T and S 
    117             CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )          
     119            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
    118120         ENDIF 
    119          tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
    120          sshn (:,:)     = sshb(:,:)    
    121          un   (:,:,:)   = ub  (:,:,:) 
    122          vn   (:,:,:)   = vb  (:,:,:) 
    123          hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level 
    124          CALL div_hor( 0 )                    ! compute interior hdivn value   
    125 !!gm                                    hdivn(:,:,:) = 0._wp 
     121         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     122         ssh (:,:,Kmm)     = ssh(:,:,Kbb)    
     123         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     124         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
     125         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level 
     126         CALL div_hor( 0, Kbb, Kmm )         ! compute interior hdiv value   
     127!!gm                                    hdiv(:,:,:) = 0._wp 
    126128 
    127129!!gm POTENTIAL BUG : 
    128 !!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
     130!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
    129131!!             as well as gdept and gdepw....   !!!!!  
    130132!!      ===>>>>   probably a call to domvvl initialisation here.... 
     
    136138!            ALLOCATE( zuvd(jpi,jpj,jpk,2) ) 
    137139!            CALL dta_uvd( nit000, zuvd ) 
    138 !            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:) 
    139 !            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:) 
     140!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb) 
     141!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb) 
    140142!            DEALLOCATE( zuvd ) 
    141143!         ENDIF 
    142144         ! 
    143145!!gm This is to be changed !!!! 
    144 !         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 
     146!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here 
    145147!         IF( .NOT.ln_linssh ) THEN 
    146148!            DO jk = 1, jpk 
    147 !               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
     149!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 
    148150!            END DO 
    149151!         ENDIF 
     
    155157      ! Do it whatever the free surface method, these arrays being eventually used 
    156158      ! 
    157       un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp 
    158       ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
     159      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp 
     160      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp 
    159161      ! 
    160 !!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked 
     162!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked 
    161163      DO jk = 1, jpkm1 
    162164         DO jj = 1, jpj 
    163165            DO ji = 1, jpi 
    164                un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 
    165                vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     166               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
     167               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
    166168               ! 
    167                ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 
    168                vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 
     169               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 
     170               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 
    169171            END DO 
    170172         END DO 
    171173      END DO 
    172174      ! 
    173       un_b(:,:) = un_b(:,:) * r1_hu_n(:,:) 
    174       vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:) 
     175      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 
     176      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm) 
    175177      ! 
    176       ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 
    177       vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 
     178      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb) 
     179      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb) 
    178180      ! 
    179181   END SUBROUTINE istate_init 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/sbcfwb.F90

    r11931 r12077  
    1717   USE dom_oce        ! ocean space and time domain 
    1818   USE sbc_oce        ! surface ocean boundary condition 
     19   USE isf_oce        ! ice shelf melting contribution 
    1920   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass 
    2021   USE phycst         ! physical constants 
    2122   USE sbcrnf         ! ocean runoffs 
    22    USE isf            ! ice shelf melting contribution 
    2323   USE sbcssr         ! Sea-Surface damping terms 
    24    USE tradmp         ! 
    2524   ! 
    2625   USE in_out_manager ! I/O manager 
     
    4948CONTAINS 
    5049 
    51    SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc ) 
     50   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc, Kmm ) 
    5251      !!--------------------------------------------------------------------- 
    5352      !!                  ***  ROUTINE sbc_fwb  *** 
     
    6665      INTEGER, INTENT( in ) ::   kn_fsbc  !  
    6766      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index 
     67      INTEGER, INTENT( in ) ::   Kmm      ! ocean time level index 
    6868      ! 
    6969      INTEGER  ::   inum, ikty, iyear     ! local integers 
     
    157157            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow 
    158158                                                      ! sum over the global domain 
    159             a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) ) 
     159            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rau0 ) ) 
    160160            a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s 
    161161!!gm        !                                                      !!bug 365d year  
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/tradmp.F90

    r11889 r12077  
    7272 
    7373 
    74    SUBROUTINE tra_dmp( kt ) 
     74   SUBROUTINE tra_dmp( kt, Kbb, Kmm, pts, Krhs ) 
    7575      !!---------------------------------------------------------------------- 
    7676      !!                   ***  ROUTINE tra_dmp  *** 
     
    9090      !! ** Action  : - tsa: tracer trends updated with the damping trend 
    9191      !!---------------------------------------------------------------------- 
    92       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     92      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step index 
     93      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices 
     94      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation 
    9395      ! 
    9496      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
     
    101103      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    102104         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )  
    103          ztrdts(:,:,:,:) = tsa(:,:,:,:)  
     105         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs)  
    104106      ENDIF 
    105107      !                           !==  input T-S data at kt  ==! 
     
    113115               DO jj = 2, jpjm1 
    114116                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    115                      tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) ) 
     117                     pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           & 
     118                        &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) ) 
    116119                  END DO 
    117120               END DO 
     
    124127               DO ji = fs_2, fs_jpim1   ! vector opt. 
    125128                  IF( avt(ji,jj,jk) <= avt_c ) THEN 
    126                      tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    127                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    128                      tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
    129                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     129                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
     130                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
     131                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
     132                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
    130133                  ENDIF 
    131134               END DO 
     
    137140            DO jj = 2, jpjm1 
    138141               DO ji = fs_2, fs_jpim1   ! vector opt. 
    139                   IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN 
    140                      tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    141                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    142                      tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
    143                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     142                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN 
     143                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
     144                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
     145                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
     146                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
    144147                  ENDIF 
    145148               END DO 
     
    150153      ! 
    151154      IF( l_trdtra )   THEN       ! trend diagnostic 
    152          ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:) 
    153          CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) ) 
    154          CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) ) 
     155         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:) 
     156         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) ) 
     157         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) ) 
    155158         DEALLOCATE( ztrdts )  
    156159      ENDIF 
    157160      !                           ! Control print 
    158       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' dmp  - Ta: ', mask1=tmask,   & 
    159          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     161      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   & 
     162         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    160163      ! 
    161164      IF( ln_timing )   CALL timing_stop('tra_dmp') 
     
    179182      REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation 
    180183      READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 
    181 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 
     184901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' ) 
    182185      ! 
    183186      REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation 
    184187      READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 
    185 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 
     188902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' ) 
    186189      IF(lwm) WRITE ( numond, namtra_dmp ) 
    187190      ! 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/usrdef_sbc.F90

    r11889 r12077  
    4141CONTAINS 
    4242 
    43    SUBROUTINE usrdef_sbc_oce( kt ) 
     43   SUBROUTINE usrdef_sbc_oce( kt, Kbb ) 
    4444      !!--------------------------------------------------------------------- 
    4545      !!                    ***  ROUTINE usr_def_sbc  *** 
     
    5656      !!---------------------------------------------------------------------- 
    5757      INTEGER, INTENT(in) ::   kt   ! ocean time step 
     58      INTEGER, INTENT(in) ::   Kbb  ! ocean time index 
    5859      !!--------------------------------------------------------------------- 
    5960      ! 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP/EXPREF/namelist_cfg

    r12068 r12077  
    150150   ! ---------------- ice shelf load ------------------------------- 
    151151   ! 
    152    cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc) 
    153152   ! 
    154153   ! ---------------- ice shelf melt formulation ------------------------------- 
  • NEMO/branches/2019/UKMO_MERGE_2019/tests/demo_cfgs.txt

    r10516 r12077  
    11CANAL OCE 
    22ISOMIP OCE 
     3ISOMIP+ OCE 
    34LOCK_EXCHANGE OCE 
    45OVERFLOW OCE 
Note: See TracChangeset for help on using the changeset viewer.