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

Changeset 13080


Ignore:
Timestamp:
2020-06-09T18:45:11+02:00 (4 years ago)
Author:
dancopsey
Message:

Merge in existing fix_cpl branch.

Location:
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt
Files:
25 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/cfgs/SHARED/field_def_nemo-ice.xml

    r11575 r13080  
    158158          <field id="hfxcndtop"    long_name="Net conductive heat flux at the ice surface (neg = ice cooling)"  standard_name="conductive_heat_flux_at_sea_ice_surface"    unit="W/m2" />  
    159159          <field id="hfxcndbot"    long_name="Net conductive heat flux at the ice bottom (neg = ice cooling)"   standard_name="conductive_heat_flux_at_sea_ice_bottom"     unit="W/m2" /> 
     160          <field id="hfxcndcpl"    long_name="Conductive heat flux coming through the coupler"                  standard_name="conductive_heat_flux_from_coupler"          unit="W/m2" /> 
    160161 
    161162          <!-- diags --> 
     
    245246     <field id="iceconc_cat"  long_name="Sea-ice concentration per category"                unit=""        /> 
    246247          <field id="icethic_cat"  long_name="Sea-ice thickness per category"                    unit="m"       detect_missing_value="true" /> 
     248          <field id="icevol_cat"   long_name="Sea-ice volume per category"                       unit="m"       detect_missing_value="true" /> 
    247249          <field id="snwthic_cat"  long_name="Snow thickness per category"                       unit="m"       detect_missing_value="true" /> 
    248250          <field id="icesalt_cat"  long_name="Sea-Ice Bulk salinity per category"                unit="g/kg"    detect_missing_value="true" /> 
     
    253255          <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" />  
    254256          <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        />  
     257          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        />  
     258          <field id="icelhpnd_cat" long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" />  
    255259          <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        /> 
    256260          <field id="iceage_cat"   long_name="Ice age per category"                              unit="days"    detect_missing_value="true" /> 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/cfgs/SHARED/namelist_ice_ref

    r11649 r13080  
    177177   ln_pnd           = .false.         !  activate melt ponds or not 
    178178     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
     179       rn_pnd_min   = 0.15            !  Minimum ice fraction that contributes to melt ponds 
     180       rn_pnd_max   = 0.7             !  Maximum ice fraction that contributes to melt ponds 
     181       ln_pnd_totfrac = .false.       !  Use total ice fraction instead of category ice fraction in melt pond contribution fracton 
     182       ln_pnd_overflow = .false.      !  Allow ponds to overflow and have vertical flushing 
     183       ln_pnd_lids  = .false.         !  Melt ponds can have frozen lids 
     184       ln_use_pndmass = .true.        !  Use melt pond mass flux diagnostic, passing it to the ocean for emp 
    179185     ln_pnd_CST     = .false.         !  activate constant  melt ponds 
    180186       rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
     
    186192!------------------------------------------------------------------------------ 
    187193   ln_iceini        = .true.          !  activate ice initialization (T) or not (F) 
    188    ln_iceini_file   = .false.         !  netcdf file provided for initialization (T) or not (F) 
     194   nn_iceini_file   = 0               !        0 = Initialise sea ice based on SSTs 
     195                                      !        1 = Initialise sea ice from single category netcdf file 
     196                                      !        2 = Initialise sea ice from multi category restart file 
    189197   rn_thres_sst     =   2.0           !  max temp. above Tfreeze with initial ice = (sst - tfreeze) 
    190198   rn_hti_ini_n     =   3.0           !  initial ice thickness       (m), North 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/doc/namelists/namthd_pnd

    r11536 r13080  
    44   ln_pnd           = .false.         !  activate melt ponds or not 
    55     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
     6       rn_pnd_min   = 0.15            !  Minimum ice fraction that contributes to melt ponds 
     7       rn_pnd_max   = 0.7             !  Maximum ice fraction that contributes to melt ponds 
     8       ln_pnd_totfrac = .false.       !  Use total ice fraction instead of category ice fraction in melt pond contribution fracton 
     9       ln_pnd_overflow = .false.      !  Allow ponds to overflow and have vertical flushing 
     10       ln_pnd_lids  = .false.         !  Melt ponds can have frozen lids 
     11       ln_use_pndmass = .true.        !  Use melt pond mass flux diagnostic, passing it to the ocean for emp 
    612     ln_pnd_CST     = .false.         !  activate constant  melt ponds 
    713       rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice.F90

    r11715 r13080  
    7070   !! a_ip        |      -      |    Ice pond concentration       |       | 
    7171   !! v_ip        |      -      |    Ice pond volume per unit area| m     | 
     72   !! lh_ip       !    lh_ip_1d !    Ice pond lid thickness       ! m     ! 
    7273   !!                                                                     | 
    7374   !!-------------|-------------|---------------------------------|-------| 
     
    195196   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1) 
    196197   LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo 
     198   REAL(wp), PUBLIC ::   rn_pnd_min       !: Minimum ice fraction that contributes to melt ponds 
     199   REAL(wp), PUBLIC ::   rn_pnd_max       !: Maximum ice fraction that contributes to melt ponds 
     200   LOGICAL,  PUBLIC ::   ln_pnd_totfrac   !: Use total ice fraction instead of category ice fraction 
     201                                          !: when determining ice fraction that contributes to melt ponds 
     202   LOGICAL,  PUBLIC ::   ln_pnd_overflow  !: Allow ponds to overflow and have vertical flushing 
     203   LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Melt ponds can have frozen lids 
     204   LOGICAL,  PUBLIC ::   ln_use_pndmass   !: Use melt pond mass flux diagnostic, passing it to the ocean for emp 
    197205 
    198206   !                                     !!** ice-diagnostics namelist (namdia) ** 
     
    295303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_i       !: Ice thickness                           (m) 
    296304   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i       !: Ice fractional areas (concentration) 
     305   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple    !: Ice fractional area at last coupling time 
    297306   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_i       !: Ice volume per unit area                (m) 
    298307   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_s       !: Snow volume per unit area               (m) 
     
    331340   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area      [m] 
    332341   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i) 
     342   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_eff   !: melt pond effective fraction (not covered up by lid) (a_ip/a_i) 
    333343   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m] 
     344   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   lh_ip      !: melt pond lid thickness                  [m] 
    334345 
    335346   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration 
     
    435446 
    436447      ii = ii + 1 
     448      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) ) 
     449 
     450      ii = ii + 1 
    437451      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                   & 
    438452         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , st_i(jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,  & 
     
    448462 
    449463      ii = ii + 1 
    450       ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 
     464      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  & 
     465         &      lh_ip(jpi,jpj,jpl), a_ip_eff(jpi,jpj,jpl) , STAT = ierr(ii) ) 
    451466 
    452467      ii = ii + 1 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice1d.F90

    r11715 r13080  
    128128   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !: 
    129129   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !: 
     130   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_eff_1d   !: 
     131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   lh_ip_1d      !: Ice pond lid thickness   [m] 
    130132 
    131133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     
    157159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d 
    158160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d  
     161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lh_ip_2d  
    159162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d  
    160163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_2d 
     
    208211         &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &     
    209212         &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  & 
    210          &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) ,                   & 
    211          &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   & 
     213         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , lh_ip_1d(jpij) ,  & 
     214         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) , a_ip_eff_1d(jpij) ,                               & 
    212215         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
    213216      ! 
     
    226229      ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  & 
    227230         &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  & 
    228          &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                      & 
     231         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl),  & 
    229232         &      STAT=ierr(ii) ) 
    230233 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icealb.F90

    r11715 r13080  
    9696      LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo 
    9797      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area) 
     98                                                                  !  This is the effective fraction not covered up by a pond lid 
    9899      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth 
    99100      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn.F90

    r11715 r13080  
    101101      ELSEWHERE 
    102102         h_ip(:,:,:) = 0._wp 
     103         lh_ip(:,:,:) = 0._wp 
    103104      END WHERE 
    104105      ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_adv.F90

    r11715 r13080  
    8484         !                             !-----------------------! 
    8585         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 
    86             &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     86            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i ) 
    8787         !                             !-----------------------! 
    8888      CASE( np_advPRA )                ! PRATHER scheme        ! 
    8989         !                             !-----------------------! 
    9090         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, & 
    91             &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     91            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i ) 
    9292      END SELECT 
    9393 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_adv_pra.F90

    r11715 r13080  
    4444   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction 
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxlhp, sylhp, sxxlhp, syylhp, sxylhp   ! melt pond lid thickness 
    4647 
    4748   !! * Substitutions 
     
    5556 
    5657   SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  & 
    57       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    5859      !!---------------------------------------------------------------------- 
    5960      !!                **  routine ice_dyn_adv_pra  ** 
     
    7879      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    7980      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     81      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    8082      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8183      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    8890      REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw 
    8991      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    90       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     92      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0lhp 
    9193      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    9294      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     
    129131            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    130132            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     133            z0lhp(:,:,jl)  = plh_ip(:,:,jl) * e1e2t(:,:)   ! Melt pond lid thickness 
    131134         ENDIF 
    132135      END DO 
     
    167170               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    168171               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     172               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness 
     173               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 
    169174            ENDIF 
    170175         END DO 
     
    202207               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    203208               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     209               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness 
     210               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)  
    204211            ENDIF 
    205212         END DO 
     
    225232            pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    226233            pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     234            plh_ip(:,:,jl) = z0lhp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    227235         ENDIF 
    228236      END DO 
     
    231239      !     Remove negative values (conservation is ensured) 
    232240      !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    233       CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     241      CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    234242      ! 
    235243      ! --- Ensure snow load is not too big --- ! 
     
    653661         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
    654662         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     663         &      sxlhp(jpi,jpj,jpl) , sylhp (jpi,jpj,jpl), sxxlhp (jpi,jpj,jpl), syylhp (jpi,jpj,jpl), sxylhp (jpi,jpj,jpl),   & 
    655664         ! 
    656665         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & 
     
    765774               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 
    766775               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 
     776               !                                                     ! melt pond lid thickness 
     777               CALL iom_get( numrir, jpdom_autoglo, 'sxlhp' , sxlhp  ) 
     778               CALL iom_get( numrir, jpdom_autoglo, 'sylhp' , sylhp  ) 
     779               CALL iom_get( numrir, jpdom_autoglo, 'sxxlhp', sxxlhp ) 
     780               CALL iom_get( numrir, jpdom_autoglo, 'syylhp', syylhp ) 
     781               CALL iom_get( numrir, jpdom_autoglo, 'sxylhp', sxylhp ) 
    767782            ENDIF 
    768783            ! 
     
    782797               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction 
    783798               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume 
     799               sxlhp  = 0._wp  ;   sylhp  = 0._wp  ;   sxxlhp  = 0._wp  ;   syylhp  = 0._wp  ;   sxylhp  = 0._wp  ! melt pond lid thickness 
    784800            ENDIF 
    785801         ENDIF 
     
    862878            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 
    863879            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 
     880            !                                                        ! melt pond lid thickness 
     881            CALL iom_rstput( iter, nitrst, numriw, 'sxlhp' , sxlhp  ) 
     882            CALL iom_rstput( iter, nitrst, numriw, 'sylhp' , sylhp  ) 
     883            CALL iom_rstput( iter, nitrst, numriw, 'sxxlhp', sxxlhp ) 
     884            CALL iom_rstput( iter, nitrst, numriw, 'syylhp', syylhp ) 
     885            CALL iom_rstput( iter, nitrst, numriw, 'sxylhp', sxylhp ) 
    864886         ENDIF 
    865887         ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_adv_umx.F90

    r11715 r13080  
    6060 
    6161   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    62       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     62      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    6363      !!---------------------------------------------------------------------- 
    6464      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
     
    8585      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration 
    8686      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     87      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    8788      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8889      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    334335            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
    335336               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
     337            ! lid thickness 
     338            zamsk = 0._wp 
     339            zhvar(:,:,:) = plh_ip(:,:,:) * z1_aip(:,:,:) 
     340            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     341               &                                      zhvar, plh_ip, zua_ups, zva_ups ) 
     342             
    336343         ENDIF 
    337344         ! 
     
    350357         ! Remove negative values (conservation is ensured) 
    351358         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    352          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     359         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    353360         ! 
    354361         ! Make sure ice thickness is not too big 
     
    15191526      !! 
    15201527      !! ** Purpose : Thickness correction in case advection scheme creates 
    1521       !!              abnormally tick ice or snow 
     1528      !!              abnormally thick ice or snow 
    15221529      !! 
    15231530      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_rdgrft.F90

    r11715 r13080  
    503503      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1 
    504504      REAL(wp)                  ::   airft1, oirft1, aprft1 
    505       REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges 
    506       REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice 
     505      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, lhprdg  ! area etc of new ridges 
     506      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft  ! area etc of rafted ice 
    507507      ! 
    508508      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges 
     
    578578                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
    579579                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg 
     580                  lhprdg(ji) = lh_ip_2d(ji,jl1) * afrdg 
    580581                  aprft1     = a_ip_2d(ji,jl1) * afrft 
    581582                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 
    582583                  vprft (ji) = v_ip_2d(ji,jl1) * afrft 
     584                  lhprft(ji) = lh_ip_2d(ji,jl1) * afrft 
    583585               ENDIF 
    584586 
     
    610612                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    611613                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
     614                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - lhprdg(ji) - lhprft(ji) 
    612615               ENDIF 
    613616            ENDIF 
     
    706709                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         &  
    707710                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   ) 
     711                     lh_ip_2d (ji,jl2) = lh_ip_2d(ji,jl2) + (   lhprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
     712                        &                                   + lhprft (ji) * rn_fpndrft * zswitch(ji)   ) 
    708713                  ENDIF 
    709714                   
     
    736741      !---------------- 
    737742      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    738       CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     743      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d ) 
    739744      ! 
    740745   END SUBROUTINE rdgrft_shift 
     
    848853         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    849854         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     855         CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 
    850856         DO jl = 1, jpl 
    851857            DO jk = 1, nlay_s 
     
    874880         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    875881         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     882         CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 
    876883         DO jl = 1, jpl 
    877884            DO jk = 1, nlay_s 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceistate.F90

    r11715 r13080  
    4141   !                             !! ** namelist (namini) ** 
    4242   LOGICAL, PUBLIC  ::   ln_iceini        !: Ice initialization or not 
    43    LOGICAL, PUBLIC  ::   ln_iceini_file   !: Ice initialization from 2D netcdf file 
     43   INTEGER, PUBLIC  ::   nn_iceini_file   ! Ice initialization: 
     44                                  !        0 = Initialise sea ice based on SSTs 
     45                                  !        1 = Initialise sea ice from single category netcdf file 
     46                                  !        2 = Initialise sea ice from multi category restart file 
    4447   REAL(wp) ::   rn_thres_sst 
    4548   REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n 
     
    4851   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s 
    4952   ! 
    50    !                              ! if ln_iceini_file = T 
     53   !                              ! if nn_iceini_file = 1 
    5154   INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read 
    5255   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m) 
     
    155158      a_ip     (:,:,:) = 0._wp 
    156159      v_ip     (:,:,:) = 0._wp 
     160      lh_ip    (:,:,:) = 0._wp 
    157161      a_ip_frac(:,:,:) = 0._wp 
     162      a_ip_eff (:,:,:) = 0._wp 
    158163      h_ip     (:,:,:) = 0._wp 
    159164      ! 
     
    167172      IF( ln_iceini ) THEN 
    168173         !                             !---------------! 
    169          IF( ln_iceini_file )THEN      ! Read a file   ! 
     174         IF( nn_iceini_file == 1 )THEN      ! Read a file   ! 
    170175            !                          !---------------! 
    171176            WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp 
     
    362367            a_ip_frac(:,:,:) = 0._wp 
    363368         END WHERE 
     369         a_ip_eff(:,:,:) = a_ip_frac(:,:,:) 
    364370         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    365371           
     
    466472      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read 
    467473      ! 
    468       NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, & 
     474      NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, & 
    469475         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, & 
    470476         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & 
     
    493499         WRITE(numout,*) '   Namelist namini:' 
    494500         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini 
    495          WRITE(numout,*) '      ice initialization from a netcdf file            ln_iceini_file = ', ln_iceini_file 
     501         WRITE(numout,*) '      ice initialization from a netcdf file            nn_iceini_file = ', nn_iceini_file 
    496502         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst 
    497          IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN 
     503         IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN 
    498504            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s  
    499505            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s 
     
    508514      ENDIF 
    509515      ! 
    510       IF( ln_iceini_file ) THEN                      ! Ice initialization using input file 
     516      IF( nn_iceini_file == 1 ) THEN                      ! Ice initialization using input file 
    511517         ! 
    512518         ! set si structure 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceitd.F90

    r11715 r13080  
    411411      CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 
    412412      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
     413      CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 
    413414      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
    414415      DO jl = 1, jpl 
     
    483484                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    484485                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
     486                  !                                               
     487                  ztrans          = lh_ip_2d(ji,jl1) * zworka(ji)     ! Pond lid thickness 
     488                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - ztrans 
     489                  lh_ip_2d(ji,jl2) = lh_ip_2d(ji,jl2) + ztrans 
    485490               ENDIF 
    486491               ! 
     
    527532      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
    528533      !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
    529       CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     534      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d ) 
    530535 
    531536      ! at_i must be <= rn_amax 
     
    555560      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 
    556561      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
     562      CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 
    557563      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
    558564      DO jl = 1, jpl 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icerst.F90

    r11715 r13080  
    132132      CALL iom_rstput( iter, nitrst, numriw, 'a_ip' , a_ip  ) 
    133133      CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip  ) 
     134      CALL iom_rstput( iter, nitrst, numriw, 'lh_ip', lh_ip ) 
    134135      ! Snow enthalpy 
    135136      DO jk = 1, nlay_s  
     
    171172      INTEGER           ::   jk 
    172173      LOGICAL           ::   llok 
    173       INTEGER           ::   id0, id1, id2, id3, id4   ! local integer 
     174      INTEGER           ::   id0, id1, id2, id3, id4, id5   ! local integer 
    174175      CHARACTER(len=25) ::   znam 
    175176      CHARACTER(len=2)  ::   zchar, zchar1 
     
    250251            v_ip(:,:,:) = 0._wp 
    251252         ENDIF 
     253         a_ip_eff(:,:,:) = a_ip(:,:,:) 
     254         ! melt pond lids 
     255         id5 = iom_varid( numrir, 'lh_ip' , ldstop = .FALSE. ) 
     256         IF( id5 > 0 ) THEN 
     257            CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip) 
     258         ELSE 
     259            lh_ip(:,:,:) = 0._wp 
     260         ENDIF 
    252261         ! fields needed for Met Office (Jules) coupling 
    253262         IF( ln_cpl ) THEN 
     
    274283         CALL ice_istate( nit000 ) 
    275284         ! 
    276          IF( .NOT.ln_iceini .OR. .NOT.ln_iceini_file ) & 
    277             &   CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and ln_iceini_file=T') 
     285         IF( .NOT.ln_iceini .OR. nn_iceini_file == 0 ) & 
     286            &   CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and nn_iceini_file=0') 
    278287         ! 
    279288      ENDIF 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icesbc.F90

    r11715 r13080  
    132132 
    133133      ! --- cloud-sky and overcast-sky ice albedos --- ! 
    134       CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) 
     134      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os ) 
    135135 
    136136      ! albedo depends on cloud fraction because of non-linear spectral effects 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icestp.F90

    r11715 r13080  
    252252      ! 
    253253      !                                ! Initial sea-ice state 
    254       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
    255          CALL ice_istate_init 
    256          CALL ice_istate( nit000 ) 
    257       ELSE                                    ! start from a restart file 
    258          CALL ice_rst_read 
     254 
     255      CALL ice_istate_init 
     256      IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 
     257         CALL ice_rst_read                      ! start from a restart file 
     258      ELSE 
     259         CALL ice_istate( nit000 )              ! start from rest: sea-ice deduced from sst 
    259260      ENDIF 
    260261      CALL ice_var_glo2eqv 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd.F90

    r11715 r13080  
    355355         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    356356         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
     357         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 
     358         CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d     (1:npti), lh_ip     (:,:,kl) ) 
    357359         ! 
    358360         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            ) 
     
    461463         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    462464         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
     465         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 
     466         CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d    (1:npti), lh_ip    (:,:,kl) ) 
    463467         ! 
    464468         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_pnd.F90

    r11715 r13080  
    3838   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012) 
    3939 
     40   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp 
     41 
    4042   !! * Substitutions 
    4143#  include "vectopt_loop_substitute.h90" 
     
    8991         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 
    9092            a_ip_frac_1d(ji) = rn_apnd 
     93            a_ip_eff_1d(ji)  = rn_apnd 
    9194            h_ip_1d(ji)      = rn_hpnd     
    9295            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji) 
    9396         ELSE 
    9497            a_ip_frac_1d(ji) = 0._wp 
     98            a_ip_eff_1d(ji)  = 0._wp 
    9599            h_ip_1d(ji)      = 0._wp     
    96100            a_ip_1d(ji)      = 0._wp 
     
    129133      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio 
    130134      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
    131       ! 
     135      REAL(wp), PARAMETER ::   pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 
     136      REAL(wp), PARAMETER ::   pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation 
     137      ! 
     138      REAL(wp) ::   tot_mlt          ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 
    132139      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding 
    133       REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding 
     140      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change) 
    134141      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    135142      REAL(wp) ::   z1_rhow          ! inverse freshwater density 
    136143      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio 
     144      REAL(wp) ::   z1_rhoi          ! inverse ice density 
    137145      REAL(wp) ::   zfac, zdum 
    138       ! 
    139       INTEGER  ::   ji   ! loop indices 
     146      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing 
     147      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing 
     148      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable) 
     149      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m) 
     150      REAL(wp) ::   dh_ip_over       ! Pond thickness change due to leaking 
     151      REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean 
     152      REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond 
     153      REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice 
     154      REAL(wp) ::   Sbr              ! Brine salinity 
     155      REAL(wp), DIMENSION(nlay_i) ::   phi     ! liquid fraction 
     156      REAL(wp) ::   perm             ! Permeability of the sea ice 
     157      REAL(wp) ::   za_ip            ! Temporary pond fraction 
     158      REAL(wp) ::   max_h_diff_ts    ! Maximum meltpond depth change due to leaking or overflow (m per ts) 
     159      REAL(wp) ::   lfrac_pnd        ! The fraction of the meltpond exposed (not inder a frozen lid) 
     160       
     161      ! 
     162      INTEGER  ::   ji, jk   ! loop indices 
    140163      !!------------------------------------------------------------------- 
    141164      z1_rhow        = 1._wp / rhow  
    142165      z1_zpnd_aspect = 1._wp / zpnd_aspect 
    143166      z1_Tp          = 1._wp / zTp  
     167      z1_rhoi        = 1._wp / rhoi 
     168 
     169      ! Define time-independent field for use in refreezing 
     170      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow) 
    144171 
    145172      DO ji = 1, npti 
    146          !                                                        !--------------------------------! 
    147          IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  ! 
    148             !                                                     !--------------------------------! 
    149             !--- Remove ponds on thin ice 
     173 
     174         !                                                            !----------------------------------------------------! 
     175         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
     176            !                                                         !----------------------------------------------------! 
     177            !--- Remove ponds on thin ice or tiny ice fractions 
    150178            a_ip_1d(ji)      = 0._wp 
    151179            a_ip_frac_1d(ji) = 0._wp 
    152180            h_ip_1d(ji)      = 0._wp 
     181            lh_ip_1d(ji)     = 0._wp 
    153182            !                                                     !--------------------------------! 
    154183         ELSE                                                     ! Case ice thickness >= rn_himin ! 
    155184            !                                                     !--------------------------------! 
    156185            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step 
    157             ! 
    158             ! available meltwater for melt ponding [m, >0] and fraction 
    159             zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    160             zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    161             !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
     186 
     187            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 
     188            za_ip = a_ip_1d(ji) 
     189            IF ( za_ip < epsi06 ) za_ip = epsi06 
     190            ! 
     191            ! available meltwater for melt ponding [m, >0] and fraction  
     192            tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow  
     193            IF ( ln_pnd_totfrac ) THEN 
     194               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction 
     195            ELSE 
     196               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji)  ! Use category ice fraction  
     197            ENDIF  
     198            zdv_mlt = zfr_mlt * tot_mlt  
    162199            ! 
    163200            !--- Pond gowth ---! 
    164             ! v_ip should never be negative, otherwise code crashes 
    165             v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 
     201            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     202            ! 
     203            !--- Lid shrinking. ---! 
     204            IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 
    166205            ! 
    167206            ! melt pond mass flux (<0) 
    168             IF( zdv_mlt > 0._wp ) THEN 
    169                zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice 
     207            IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN 
     208               zfac = zdv_mlt * rhow * r1_rdtice 
    170209               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    171210               ! 
     
    177216            ! 
    178217            !--- Pond contraction (due to refreezing) ---! 
    179             v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 
     218            IF ( ln_pnd_lids ) THEN 
     219 
     220               ! Code to use if using melt pond lids 
     221               IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 
     222                  t_grad = (zTp+rt0) - t_su_1d(ji) 
     223 
     224                  ! The following equation is a rearranged form of: 
     225                  ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 
     226                  ! where: lid_thickness_start = lh_ip_1d(ji) 
     227                  !        lid_thickness_end = lh_ip_end 
     228                  !        omega_dt is a bunch of terms in the equation that do not change 
     229                  ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 
     230                  ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density).  
     231 
     232                  lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
     233                  zdh_frz = lh_ip_end - lh_ip_1d(ji) 
     234 
     235                  ! Pond shrinking 
     236                  v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 
     237 
     238                  ! Lid growing 
     239                  IF ( ln_pnd_lids )  lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 
     240               ELSE 
     241                  zdh_frz = 0._wp 
     242               END IF 
     243 
     244            ELSE 
     245               ! Code to use if not using melt pond lids 
     246               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 
     247            ENDIF 
     248 
     249            ! 
     250            ! Make sure pond volume or lid thickness has not gone negative 
     251            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp  
     252            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 
    180253            ! 
    181254            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 
     
    184257            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    185258            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji) 
     259             
     260            !--- Pond overflow and vertical flushing ---! 
     261            IF ( ln_pnd_overflow ) THEN 
     262 
     263               ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
     264               IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN 
     265                   h_ip_1d(ji) = zpnd_aspect * zfr_mlt 
     266                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     267                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     268                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     269               ENDIF 
     270 
     271               ! If pond depth exceeds half the ice thickness then reduce the pond volume 
     272               IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 
     273                   h_ip_1d(ji) = 0.5_wp * h_i_1d(ji) 
     274                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     275                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     276                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     277               ENDIF 
     278 
     279               !-- Vertical flushing of pond water --! 
     280               ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth. 
     281               ! This assumes the pond is sitting on top of the ice. 
     282               h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 
     283 
     284               ! The depth through which water percolates is the distance under the melt pond 
     285               h_percolation = h_i_1d(ji) - h_ip_1d(ji) 
     286 
     287               ! Calculate the permeability of the ice (Assur 1958) 
     288               DO jk = 1, nlay_i 
     289                   Sbr = - 1.2_wp                         & 
     290                         - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    & 
     291                         - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 & 
     292                         - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3 
     293                   phi(jk) = sz_i_1d(ji,jk)/Sbr 
     294               END DO 
     295               perm = 3.0e-08_wp * (minval(phi))**3 
     296 
     297               ! Do the drainage using Darcy's law 
     298               IF ( perm > 0._wp ) THEN 
     299                   dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number 
     300                   dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative 
     301 
     302                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
     303                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     304                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     305                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     306               ENDIF 
     307            ENDIF 
     308 
     309            ! If lid thickness is ten times greater than pond thickness then remove pond  
     310            IF ( ln_pnd_lids ) THEN            
     311               IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     312                   a_ip_1d(ji)      = 0._wp 
     313                   a_ip_frac_1d(ji) = 0._wp 
     314                   h_ip_1d(ji)      = 0._wp 
     315                   lh_ip_1d(ji)     = 0._wp 
     316                   v_ip_1d(ji)      = 0._wp 
     317               ENDIF 
     318            ENDIF 
     319 
     320            ! If any of the previous changes has removed all the ice thickness then remove ice area. 
     321            IF ( h_i_1d(ji) == 0._wp ) THEN 
     322                a_i_1d(ji) = 0._wp 
     323                h_s_1d(ji) = 0._wp 
     324            ENDIF 
     325 
    186326            ! 
    187327         ENDIF 
     328 
     329         ! Calculate how much meltpond is exposed to the atmosphere (not under a frozen lid)         
     330         IF ( lh_ip_1d(ji) < pnd_lid_min ) THEN       ! Pond lid is very thin. Expose all the pond. 
     331            lfrac_pnd = 1.0 
     332         ELSE 
     333            IF ( lh_ip_1d(ji) > pnd_lid_max ) THEN    ! Pond lid is very thick. Cover all the pond up with ice and nosw. 
     334              lfrac_pnd = 0.0 
     335            ELSE                                      ! Pond lid is of intermediate thickness. Expose part of the melt pond. 
     336              lfrac_pnd = ( lh_ip_1d(ji) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 
     337            ENDIF 
     338         ENDIF 
     339          
     340         ! Calculate the melt pond effective area 
     341         a_ip_eff_1d(ji) = a_ip_frac_1d(ji) * lfrac_pnd 
     342 
    188343      END DO 
    189344      ! 
     
    205360      INTEGER  ::   ios, ioptio   ! Local integer 
    206361      !! 
    207       NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
     362      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb,          & 
     363                            rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac,  & 
     364                            ln_use_pndmass 
    208365      !!------------------------------------------------------------------- 
    209366      ! 
     
    223380         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd 
    224381         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 
     382         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds       rn_pnd_min = ', rn_pnd_min 
     383         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds       rn_pnd_max = ', rn_pnd_max 
     384         WRITE(numout,*) '            Use total ice fraction instead of category ice fraction   ln_pnd_totfrac  = ',ln_pnd_totfrac 
     385         WRITE(numout,*) '            Allow ponds to overflow and have vertical flushing        ln_pnd_overflow = ',ln_pnd_overflow 
     386         WRITE(numout,*) '            Melt ponds can have frozen lids                           ln_pnd_lids     = ',ln_pnd_lids 
     387         WRITE(numout,*) '            Use melt pond mass flux diagnostic, passing to ocean      ln_use_pndmass  = ',ln_use_pndmass 
    225388         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST 
    226389         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_zdf_bl99.F90

    r11715 r13080  
    3131   !!---------------------------------------------------------------------- 
    3232   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    33    !! $Id$ 
     33   !! $Id: icethd_zdf_bl99.F90 10926 2019-05-03 12:32:10Z clem $ 
    3434   !! Software governed by the CeCILL license (see ./LICENSE) 
    3535   !!---------------------------------------------------------------------- 
     
    769769      ! 
    770770      ! --- calculate conduction fluxes (positive downward) 
    771  
     771      !     bottom ice conduction flux 
    772772      DO ji = 1, npti 
    773          !                                ! surface ice conduction flux 
    774          qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
    775             &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
    776          !                                ! bottom ice conduction flux 
    777          qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
     773         qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
    778774      END DO 
    779        
     775      !     surface ice conduction flux 
     776      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN 
     777         ! 
     778         DO ji = 1, npti 
     779            qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
     780               &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
     781         END DO 
     782         ! 
     783      ELSEIF( k_cnd == np_cnd_ON ) THEN 
     784         ! 
     785         DO ji = 1, npti 
     786            qcn_ice_top_1d(ji) = qcn_ice_1d(ji) 
     787         END DO 
     788         ! 
     789      ENDIF 
     790      !     surface ice temperature 
     791      IF( k_cnd == np_cnd_ON .AND. ln_cndemulate ) THEN 
     792         ! 
     793         DO ji = 1, npti 
     794            t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature 
     795               &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 
     796               &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) & 
     797               &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 
     798            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su 
     799         END DO 
     800         ! 
     801      ENDIF 
    780802      ! 
    781803      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 
     
    785807         DO ji = 1, npti 
    786808            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)  
    787          END DO 
    788          ! 
    789       ELSEIF( k_cnd == np_cnd_ON ) THEN 
    790          ! 
    791          DO ji = 1, npti 
    792             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)  
    793809         END DO 
    794810         ! 
     
    856872         t_i_1d    (1:npti,:) = ztiold        (1:npti,:) 
    857873         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti) 
    858  
    859          !!clem 
    860          ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input 
    861          !clem 
    862874      ENDIF 
    863875      ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceupdate.F90

    r11715 r13080  
    185185      ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 
    186186      !------------------------------------------------------------------ 
    187       CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
     187      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    188188      ! 
    189189      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     
    279279      IF( iom_use('hfxcndbot'  ) )   CALL iom_put( 'hfxcndbot'  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux 
    280280      IF( iom_use('hfxcndtop'  ) )   CALL iom_put( 'hfxcndtop'  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux 
     281      IF( iom_use('hfxcndcpl'  ) )   CALL iom_put( "hfxcndcpl"  , SUM( qcn_ice * a_i_b, dim=3 ) )       ! Conduction flux we are giving it 
    281282 
    282283      ! controls 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icevar.F90

    r11715 r13080  
    533533               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
    534534               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     535               lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj) 
    535536               ! 
    536537            END DO 
     
    555556 
    556557 
    557    SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     558   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    558559      !!------------------------------------------------------------------- 
    559560      !!                   ***  ROUTINE ice_var_zapneg *** 
     
    570571      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    571572      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     573      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    572574      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    573575      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    637639      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    638640      !                                                        but it does not change conservation, so keep it this way is ok 
     641      WHERE( plh_ip (:,:,:) < 0._wp )   plh_ip (:,:,:) = 0._wp 
    639642      ! 
    640643   END SUBROUTINE ice_var_zapneg 
    641644 
    642645 
    643    SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     646   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    644647      !!------------------------------------------------------------------- 
    645648      !!                   ***  ROUTINE ice_var_roundoff *** 
     
    654657      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    655658      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     659      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    656660      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    657661      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    668672         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    669673         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     674         WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 )    plh_ip(1:npti,:)   = 0._wp   ! lh_ip must be >= 0 
    670675      ENDIF 
    671676      ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icewri.F90

    r11715 r13080  
    149149 
    150150      ! --- category-dependent fields --- ! 
     151      IF( iom_use('icelhpnd_cat') )   CALL iom_put( 'icelhpnd_cat', lh_ip          * zmsk00l                                   ) ! melt pond lid thickness for categories 
    151152      IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0% 
    152153      IF( iom_use('iceconc_cat' ) )   CALL iom_put( 'iceconc_cat' , a_i            * zmsk00l                                   ) ! area for categories 
     
    165166      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories 
    166167      IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories 
     168 
     169      IF( iom_use('icevol_cat'  ) )   CALL iom_put( "icevol_cat"  , v_i            * zmsk00l                                   ) ! volume for categories 
     170      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories 
    167171 
    168172      !------------------ 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/OCE/BDY/bdyice.F90

    r11715 r13080  
    9898                 &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                & 
    9999                 &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      ) 
     100            CALL lbc_lnk_multi( 'bdyice', lh_ip,'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      ) 
    100101            ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk 
    101102            CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
     
    163164            a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration 
    164165            h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth 
     166            lh_ip(ji,jj, jl) = ( lh_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond lid thickness 
    165167            ! 
    166168            sz_i(ji,jj,:,jl) = s_i(ji,jj,jl) 
     
    170172               a_ip(ji,jj,jl) = 0._wp 
    171173               h_ip(ji,jj,jl) = 0._wp 
     174               lh_ip(ji,jj,jl) = 0._wp 
    172175            ENDIF 
    173176            ! 
     
    231234               a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) 
    232235               h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl) 
     236               lh_ip(ji,jj, jl) = lh_ip(ib,jb, jl) 
    233237               ! 
    234238               sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) 
     
    280284               a_ip(ji,jj,  jl) = 0._wp 
    281285               v_ip(ji,jj,  jl) = 0._wp 
     286               lh_ip(ji,jj, jl) = 0._wp 
    282287               t_su(ji,jj,  jl) = rt0 
    283288               t_s (ji,jj,:,jl) = rt0 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/OCE/SBC/sbc_ice.F90

    r11715 r13080  
    7070   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndm_ice       !: wind speed module at T-point                 [m/s] 
    7171   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: wind speed module at T-point                 [m/s] 
    72    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tsfc_ice       !: sea ice surface skin temperature (on categories) 
    7372#endif 
    7473 
     
    9392    
    9493   ! already defined in ice.F90 for SI3 
    95    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i 
     94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i               ! Sea ice fraction on categories 
     95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i_last_couple   ! Sea ice fraction on categories at the last coupling point 
    9696   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  h_i, h_s 
    9797 
     
    132132         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce   (jpi,jpj)     ,   & 
    133133         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce    (jpi,jpj)     ,   & 
    134          &      emp_ice (jpi,jpj)     , tsfc_ice (jpi,jpj,jpl) , sstfrz     (jpi,jpj)     , STAT= ierr(2) ) 
     134         &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , STAT= ierr(2) ) 
    135135#endif 
    136136 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/OCE/SBC/sbccpl.F90

    r11715 r13080  
    4848   USE lib_mpp        ! distribued memory computing library 
    4949   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     50 
     51#if defined key_oasis3  
     52   USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut  
     53#endif  
    5054 
    5155   IMPLICIT NONE 
     
    152156   INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level  
    153157   INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for semi-implicit coupling of atmos-ice fluxes) 
    154    INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area 
     158   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area fraction 
    155159   INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness 
    156160   INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity 
     
    159163 
    160164   INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent  
     165 
     166#if ! defined key_oasis3  
     167   ! Dummy variables to enable compilation when oasis3 is not being used  
     168   INTEGER                    ::   OASIS_Sent        = -1  
     169   INTEGER                    ::   OASIS_SentOut     = -1  
     170   INTEGER                    ::   OASIS_ToRest      = -1  
     171   INTEGER                    ::   OASIS_ToRestOut   = -1  
     172#endif  
    161173 
    162174   !                                  !!** namelist namsbc_cpl ** 
     
    184196   LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models 
    185197                                          !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
     198   LOGICAL     ::   ln_scale_ice_fluxes   ! Scale sea ice fluxes by the sea ice fractions at the previous coupling point 
    186199   TYPE ::   DYNARR      
    187200      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3    
     
    250263      NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   &  
    251264         &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  &  
     265         &                  ln_scale_ice_fluxes,                                                       & 
    252266         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   &  
    253267         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   &  
     
    279293      ENDIF 
    280294      IF( lwp .AND. ln_cpl ) THEN                        ! control print 
     295         WRITE(numout,*)'  ln_scale_ice_fluxes                 = ', ln_scale_ice_fluxes 
    281296         WRITE(numout,*)'  received fields (mutiple ice categogies)' 
    282297         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')' 
     
    815830      END SELECT 
    816831 
     832      ! Initialise ice fractions from last coupling time to zero 
     833       a_i_last_couple(:,:,:) = 0._wp 
     834 
     835 
    817836      !                                                      ! ------------------------- !  
    818837      !                                                      !      Ice Meltponds        !  
     
    12431262      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1) 
    12441263      ! 
    1245       !                                                      ! ================== ! 
    1246       !                                                      !   ice skin temp.   ! 
    1247       !                                                      ! ================== ! 
    1248 #if defined key_si3 
    1249       ! needed by Met Office 
    1250       IF( srcv(jpr_ts_ice)%laction ) THEN  
    1251          WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   tsfc_ice(:,:,:) = 0.0  
    1252          ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   tsfc_ice(:,:,:) = -60. 
    1253          ELSEWHERE                                        ;   tsfc_ice(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) 
    1254          END WHERE 
    1255       ENDIF  
    1256 #endif 
    12571264      !                                                      ! ========================= !  
    12581265      !                                                      ! Mean Sea Level Pressure   !   (taum)  
     
    16301637      !!                   sprecip           solid precipitation over the ocean   
    16311638      !!---------------------------------------------------------------------- 
    1632       REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1] 
    1633       !                                                !!           ! optional arguments, used only in 'mixed oce-ice' case 
    1634       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo  
    1635       REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius] 
    1636       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] 
    1637       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m] 
    1638       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m] 
     1639      REAL(wp), INTENT(in)   , DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1] 
     1640      !                                                   !!           ! optional arguments, used only in 'mixed oce-ice' case or for Met-Office coupling 
     1641      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo  
     1642      REAL(wp), INTENT(in)   , DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius] 
     1643      REAL(wp), INTENT(inout), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] => inout for Met-Office 
     1644      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m] 
     1645      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m] 
    16391646      ! 
    16401647      INTEGER  ::   ji, jj, jl   ! dummy loop index 
     
    16431650      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice 
    16441651      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
    1645       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice    !!gm , zfrqsr_tr_i 
     1652      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total 
     1653      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu 
    16461654      !!---------------------------------------------------------------------- 
    16471655      ! 
     
    16631671         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
    16641672         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
    1665          zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:) 
    16661673      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp 
    16671674         zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1) 
    16681675         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:) 
    16691676         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1) 
    1670          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:) 
     1677         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)          
    16711678      CASE( 'none'      )       ! Not available as for now: needs additional coding below when computing zevap_oce  
    16721679      !                         ! since fields received are not defined with none option 
     
    16751682 
    16761683#if defined key_si3 
     1684 
     1685      ! --- evaporation over ice (kg/m2/s) --- ! 
     1686      zevap_ice_total(:,:) = 0._wp    
     1687      IF (sn_rcv_emp%clcat == 'yes') THEN 
     1688         DO jl=1,jpl 
     1689            IF (ln_scale_ice_fluxes) THEN 
     1690               zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) * a_i_last_couple(:,:,jl) 
     1691            ELSE 
     1692               zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) 
     1693            ENDIF 
     1694            zevap_ice_total(:,:) = zevap_ice_total(:,:) + zevap_ice(:,:,jl) 
     1695         ENDDO 
     1696      ELSE 
     1697         zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1 ) 
     1698         zevap_ice_total(:,:) = zevap_ice(:,:,1) 
     1699      ENDIF 
     1700 
     1701      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN 
     1702         ! For conservative case zemp_ice has not been defined yet. Do it now. 
     1703         zemp_ice(:,:) = zevap_ice_total(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:) 
     1704      END IF 
     1705 
    16771706      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing) 
    16781707      zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw ) 
     
    16831712 
    16841713      ! --- evaporation over ocean (used later for qemp) --- ! 
    1685       zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) 
    1686  
    1687       ! --- evaporation over ice (kg/m2/s) --- ! 
    1688       DO jl=1,jpl 
    1689          IF (sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) 
    1690          ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF 
    1691       ENDDO 
     1714      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) 
     1715 
     1716       
     1717 
    16921718 
    16931719      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0 
     
    17271753         sprecip (:,:)   = zsprecip (:,:) 
    17281754         tprecip (:,:)   = ztprecip (:,:) 
    1729          evap_ice(:,:,:) = zevap_ice(:,:,:) 
     1755         IF (ln_scale_ice_fluxes) THEN 
     1756            ! Convert from grid box means to sea ice means 
     1757            WHERE( a_i(:,:,:) > 0.0_wp ) evap_ice(:,:,:) = zevap_ice(:,:,:) / a_i(:,:,:) 
     1758            WHERE( a_i(:,:,:) <= 0.0_wp ) evap_ice(:,:,:) = 0.0 
     1759         ELSE 
     1760            evap_ice(:,:,:) = zevap_ice(:,:,:) 
     1761         ENDIF 
    17301762         DO jl = 1, jpl 
    17311763            devap_ice(:,:,jl) = zdevap_ice(:,:) 
     
    17741806      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
    17751807      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
    1776       IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average) 
     1808      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , zevap_ice_total(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average) 
    17771809      IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
    1778          &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average) 
     1810         &                                                        - zevap_ice_total(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average) 
    17791811      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf 
    17801812      ! 
     
    17841816      CASE( 'oce only' )         ! the required field is directly provided 
    17851817         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) 
     1818         ! For Met Office sea ice non-solar fluxes are already delt with by JULES so setting to zero 
     1819         ! here so the only flux is the ocean only one. 
     1820         zqns_ice(:,:,:) = 0._wp  
    17861821      CASE( 'conservative' )     ! the required fields are directly provided 
    17871822         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 
     
    18471882         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus )  ! solid precip over ocean + snow melting 
    18481883      zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - rLfus )  ! solid precip over ice (qevap_ice=0 since atm. does not take it into account) 
    1849 !!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * picefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap 
    1850 !!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhos  ! solid precip over ice 
    18511884       
    18521885      ! --- total non solar flux (including evap/precip) --- ! 
     
    19001933      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting 
    19011934      IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average) 
    1902       IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) & 
    1903            &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average) 
     1935      IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:)  ) & 
     1936                                                                         * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average) 
    19041937      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average) 
    19051938      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
     
    19141947      CASE( 'oce only' ) 
    19151948         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) ) 
     1949         ! For Met Office sea ice solar fluxes are already delt with by JULES so setting to zero 
     1950         ! here so the only flux is the ocean only one. 
     1951         zqsr_ice(:,:,:) = 0._wp 
    19161952      CASE( 'conservative' ) 
    19171953         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
     
    19922028            ENDDO 
    19932029         ENDIF 
     2030      CASE( 'none'      )  
     2031         zdqns_ice(:,:,:) = 0._wp 
    19942032      END SELECT 
    19952033       
     
    20072045      !                                                      ! ========================= ! 
    20082046      CASE ('coupled') 
    2009          qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) 
    2010          qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) 
     2047         IF (ln_scale_ice_fluxes) THEN 
     2048            WHERE( a_i(:,:,:) > 0.0_wp ) qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2049            WHERE( a_i(:,:,:) <= 0.0_wp ) qml_ice(:,:,:) = 0.0_wp 
     2050            WHERE( a_i(:,:,:) > 0.0_wp ) qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2051            WHERE( a_i(:,:,:) <= 0.0_wp ) qcn_ice(:,:,:) = 0.0_wp 
     2052         ELSE 
     2053            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) 
     2054            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) 
     2055         ENDIF 
    20112056      END SELECT 
    20122057      ! 
     
    20172062         ! 
    20182063         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    2019          ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77) 
     2064         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77) 
    20202065         ! 
    2021          qtr_ice_top(:,:,:) = ztri * qsr_ice(:,:,:) 
    2022          WHERE( phs(:,:,:) >= 0.0_wp )   qtr_ice_top(:,:,:) = 0._wp            ! snow fully opaque 
    2023          WHERE( phi(:,:,:) <= 0.1_wp )   qtr_ice_top(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation 
     2066         WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     2067            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
     2068         ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
     2069            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri 
     2070         ELSEWHERE                                                         ! zero when hs>0 
     2071            zqtr_ice_top(:,:,:) = 0._wp 
     2072         END WHERE 
    20242073         !      
    20252074      ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==! 
     
    20272076         !                    ! ===> here we must receive the qtr_ice_top array from the coupler 
    20282077         !                           for now just assume zero (fully opaque ice) 
    2029          qtr_ice_top(:,:,:) = 0._wp 
     2078         zqtr_ice_top(:,:,:) = 0._wp 
     2079         ! 
     2080      ENDIF 
     2081      ! 
     2082      IF( ln_mixcpl ) THEN 
     2083         DO jl=1,jpl 
     2084            qtr_ice_top(:,:,jl) = qtr_ice_top(:,:,jl) * xcplmask(:,:,0) + zqtr_ice_top(:,:,jl) * zmsk(:,:) 
     2085         ENDDO 
     2086      ELSE 
     2087         qtr_ice_top(:,:,:) = zqtr_ice_top(:,:,:) 
     2088      ENDIF 
     2089      !                                                      ! ================== ! 
     2090      !                                                      !   ice skin temp.   ! 
     2091      !                                                      ! ================== ! 
     2092      ! needed by Met Office 
     2093      IF( srcv(jpr_ts_ice)%laction ) THEN  
     2094         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   ztsu(:,:,:) = 0.0 + rt0  
     2095         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   ztsu(:,:,:) = -60. + rt0 
     2096         ELSEWHERE                                        ;   ztsu(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) + rt0 
     2097         END WHERE 
     2098         ! 
     2099         IF( ln_mixcpl ) THEN 
     2100            DO jl=1,jpl 
     2101               pist(:,:,jl) = pist(:,:,jl) * xcplmask(:,:,0) + ztsu(:,:,jl) * zmsk(:,:) 
     2102            ENDDO 
     2103         ELSE 
     2104            pist(:,:,:) = ztsu(:,:,:) 
     2105         ENDIF 
    20302106         ! 
    20312107      ENDIF 
     
    20472123      INTEGER, INTENT(in) ::   kt 
    20482124      ! 
     2125      ! 
    20492126      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    20502127      INTEGER ::   isec, info   ! local integer 
    20512128      REAL(wp) ::   zumax, zvmax 
    20522129      REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 
    2053       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4    
     2130      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4 
     2131 
    20542132      !!---------------------------------------------------------------------- 
    20552133      ! 
     
    21932271      ENDIF 
    21942272 
     2273      ! If this coupling was successful then save ice fraction for use between coupling points.  
     2274      ! This is needed for some calculations where the ice fraction at the last coupling point  
     2275      ! is needed.  
     2276      IF( info == OASIS_Sent     .OR. info == OASIS_ToRest .OR.   &  
     2277                     & info == OASIS_SentOut  .OR. info == OASIS_ToRestOut ) THEN  
     2278         IF ( sn_snd_thick%clcat == 'yes' ) THEN  
     2279           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl) 
     2280         ENDIF 
     2281      ENDIF     
     2282 
    21952283      IF( ssnd(jps_fice1)%laction ) THEN 
    21962284         SELECT CASE( sn_snd_thick1%clcat ) 
     
    22502338      !                                                      !      Ice melt ponds       !  
    22512339      !                                                      ! ------------------------- ! 
    2252       ! needed by Met Office 
     2340      ! needed by Met Office - 1) fraction of ponded ice; 2) local/actual pond depth 
    22532341      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN  
    22542342         SELECT CASE( sn_snd_mpnd%cldes)   
     
    22562344            SELECT CASE( sn_snd_mpnd%clcat )   
    22572345            CASE( 'yes' )   
    2258                ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl) 
    2259                ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)   
     2346 
     2347               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl) 
     2348               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl) 
     2349 
    22602350            CASE( 'no' )   
    22612351               ztmp3(:,:,:) = 0.0   
    22622352               ztmp4(:,:,:) = 0.0   
    22632353               DO jl=1,jpl   
    2264                  ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip(:,:,jpl)   
    2265                  ztmp4(:,:,1) = ztmp4(:,:,1) + v_ip(:,:,jpl)  
     2354                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl) 
     2355                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
    22662356               ENDDO   
    22672357            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )   
Note: See TracChangeset for help on using the changeset viewer.