Changeset 13442


Ignore:
Timestamp:
2020-08-27T17:13:15+02:00 (5 months ago)
Author:
dancopsey
Message:

Put majority of code in to get melt pond freezing to feed back into ice thickness.

Location:
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice.F90

    r13080 r13442  
    7070   !! a_ip        |      -      |    Ice pond concentration       |       | 
    7171   !! v_ip        |      -      |    Ice pond volume per unit area| m     | 
     72   !! s_ip        !    s_ip_1d  !    Ice pond salinity            ! pss   ! 
    7273   !! lh_ip       !    lh_ip_1d !    Ice pond lid thickness       ! m     ! 
    7374   !!                                                                     | 
     
    271272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: salt flux due to correction on ice thick. (residual) [pss.kg.m-2.s-1 => g.m-2.s-1] 
    272273   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
     274   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_pnd     !: salt flux due to melt pond-ocean exchange            [pss.kg.m-2.s-1 => g.m-2.s-1] 
    273275 
    274276   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth           [W.m-2] 
     
    343345   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m] 
    344346   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   lh_ip      !: melt pond lid thickness                  [m] 
     347   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   s_ip       !: melt pond salinity                       [pss] 
    345348 
    346349   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration 
     
    463466      ii = ii + 1 
    464467      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) ) 
     468         &      lh_ip(jpi,jpj,jpl), a_ip_eff(jpi,jpj,jpl), s_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 
    466469 
    467470      ii = ii + 1 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice1d.F90

    r13080 r13442  
    9292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_lam_1d 
    9393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_dyn_1d 
     94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_pnd_1d 
    9495 
    9596   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d 
     
    115116   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bog      !: Ice bottom accretion  [m] 
    116117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sub      !: Ice surface sublimation [m] 
    117    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_mlt      !: Snow melt [m] 
     118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_to_pnd   !: Ice surface melt to ponds [m] 
     119   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_mlt      !: Total snow melt [m] 
     120   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_to_pnd   !: Snow melt going to ponds [m] 
    118121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice] 
    119122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_1d        !: Ice bulk salinity [ppt] 
    120123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom 
     124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_pndfrz    !: Salinity of new ice after pond freezing 
    121125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_i_1d        !: 
    122126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_s_1d        !: 
     
    126130   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !: 
    127131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d       !: 
     132   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_ip_1d       !: Ice pond salinity 
    128133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !: 
    129134   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !: 
    130135   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_eff_1d   !: 
    131136   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   lh_ip_1d      !: Ice pond lid thickness   [m] 
     137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frz_ip        !: Heat flux from freezing melt ponds 
    132138 
    133139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     
    159165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d 
    160166   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d  
     167   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   s_ip_2d  
    161168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lh_ip_2d  
    162169   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d  
     
    210217         &      h_i_1d  (jpij) , h_ib_1d (jpij) , h_s_1d    (jpij) ,                                    &     
    211218         &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &     
    212          &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  & 
     219         &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  s_i_pndfrz(jpij) ,& 
    213220         &      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) ,                              & 
     221         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) , a_ip_eff_1d(jpij), s_ip_1d(jpij) , frz_ip(jpij) , & 
    215222         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
    216223      ! 
     
    230237         &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  & 
    231238         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl),  & 
     239         &      s_ip_2d(jpij,jpl) ,                                                              & 
    232240         &      STAT=ierr(ii) ) 
    233241 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icectl.F90

    r11715 r13442  
    9494         ! salt flux 
    9595         pdiag_fs = glob_sum( 'icectl',  & 
    96             &                         ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + & 
     96            &                         ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_pnd + & 
    9797            &                           sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) 
    9898         ! heat flux 
     
    112112         ! -- salt diag -- ! 
    113113         zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice  & 
    114             &         + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni +           & 
     114            &         + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_pnd + & 
    115115            &                                 sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) & 
    116116            &         - pdiag_fs 
     
    242242            &       wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr 
    243243         ! salt flux 
    244          pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam  
     244         pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam + sfx_pnd 
    245245         ! heat flux 
    246246         pdiag_ft =   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  &  
     
    258258         ! -- salt diag -- ! 
    259259         zdiag_salt =   ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice                                                  & 
    260             &         + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) & 
     260            &         + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam + sfx_pnd ) & 
    261261            &         - pdiag_fs 
    262262         IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel )   ll_stop_s = .TRUE. 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_adv.F90

    r13080 r13442  
    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, lh_ip, e_s, e_i ) 
     86            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, sv_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, lh_ip, e_s, e_i ) 
     91            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, sv_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

    r13080 r13442  
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume 
    4646   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxlhp, sylhp, sxxlhp, syylhp, sxylhp   ! melt pond lid thickness 
     47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsalp, sysalp, sxxsalp, syysalp, sxysalp   ! melt pond salinity 
    4748 
    4849   !! * Substitutions 
     
    5657 
    5758   SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  & 
    58       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
     59      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 
    5960      !!---------------------------------------------------------------------- 
    6061      !!                **  routine ice_dyn_adv_pra  ** 
     
    8081      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
    8182      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
     83      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_ip     ! melt pond salt content 
    8284      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8385      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    132134            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
    133135            z0lhp(:,:,jl)  = plh_ip(:,:,jl) * e1e2t(:,:)   ! Melt pond lid thickness 
     136            z0smp(:,:,jl)  = psv_ip(:,:,jl) * e1e2t(:,:)   ! Melt pond salt content 
    134137         ENDIF 
    135138      END DO 
     
    172175               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness 
    173176               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 
     177               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp)    !--- melt pond salinity 
     178               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp) 
    174179            ENDIF 
    175180         END DO 
     
    209214               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness 
    210215               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)  
     216               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp)    !--- melt pond salinity 
     217               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp) 
    211218            ENDIF 
    212219         END DO 
     
    233240            pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    234241            plh_ip(:,:,jl) = z0lhp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     242            psv_ip(:,:,jl) = z0smp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    235243         ENDIF 
    236244      END DO 
     
    239247      !     Remove negative values (conservation is ensured) 
    240248      !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    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 ) 
     249      CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 
    242250      ! 
    243251      ! --- Ensure snow load is not too big --- ! 
     
    662670         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
    663671         &      sxlhp(jpi,jpj,jpl) , sylhp (jpi,jpj,jpl), sxxlhp (jpi,jpj,jpl), syylhp (jpi,jpj,jpl), sxylhp (jpi,jpj,jpl),   & 
     672         &      sxsalp(jpi,jpj,jpl) , sysalp (jpi,jpj,jpl), sxxsalp (jpi,jpj,jpl), syysalp (jpi,jpj,jpl), sxysalp (jpi,jpj,jpl),   & 
    664673         ! 
    665674         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & 
     
    780789               CALL iom_get( numrir, jpdom_autoglo, 'syylhp', syylhp ) 
    781790               CALL iom_get( numrir, jpdom_autoglo, 'sxylhp', sxylhp ) 
     791               !                                                     ! melt pond salinity 
     792               CALL iom_get( numrir, jpdom_autoglo, 'sxsalp' , sxsalp  ) 
     793               CALL iom_get( numrir, jpdom_autoglo, 'sysalp' , sysalp  ) 
     794               CALL iom_get( numrir, jpdom_autoglo, 'sxxsalp', sxxsalp ) 
     795               CALL iom_get( numrir, jpdom_autoglo, 'syysalp', syysalp ) 
     796               CALL iom_get( numrir, jpdom_autoglo, 'sxysalp', sxysalp ) 
    782797            ENDIF 
    783798            ! 
     
    798813               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume 
    799814               sxlhp  = 0._wp  ;   sylhp  = 0._wp  ;   sxxlhp  = 0._wp  ;   syylhp  = 0._wp  ;   sxylhp  = 0._wp  ! melt pond lid thickness 
     815               sxsalp  = 0._wp  ;   sysalp  = 0._wp  ;   sxxsalp  = 0._wp  ;   syysalp  = 0._wp  ;   sxysalp  = 0._wp  ! melt pond salinity 
    800816            ENDIF 
    801817         ENDIF 
     
    884900            CALL iom_rstput( iter, nitrst, numriw, 'syylhp', syylhp ) 
    885901            CALL iom_rstput( iter, nitrst, numriw, 'sxylhp', sxylhp ) 
     902            !                                                        ! melt pond salinity 
     903            CALL iom_rstput( iter, nitrst, numriw, 'sxsalp' , sxsalp  ) 
     904            CALL iom_rstput( iter, nitrst, numriw, 'sysalp' , sysalp  ) 
     905            CALL iom_rstput( iter, nitrst, numriw, 'sxxsalp', sxxsalp ) 
     906            CALL iom_rstput( iter, nitrst, numriw, 'syysalp', syysalp ) 
     907            CALL iom_rstput( iter, nitrst, numriw, 'sxysalp', sxysalp ) 
    886908         ENDIF 
    887909         ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_adv_umx.F90

    r13080 r13442  
    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, plh_ip, pe_s, pe_i ) 
     62      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 
    6363      !!---------------------------------------------------------------------- 
    6464      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
     
    8686      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
    8787      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
     88      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_ip     ! melt pond salt content 
    8889      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8990      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    340341            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
    341342               &                                      zhvar, plh_ip, zua_ups, zva_ups ) 
     343            ! salt content 
     344            zamsk = 0._wp 
     345            zhvar(:,:,:) = psv_ip(:,:,:) * z1_aip(:,:,:) 
     346            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     347               &                                      zhvar, psv_ip, zua_ups, zva_ups ) 
    342348             
    343349         ENDIF 
     
    357363         ! Remove negative values (conservation is ensured) 
    358364         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    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 ) 
     365         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 
    360366         ! 
    361367         ! Make sure ice thickness is not too big 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_rdgrft.F90

    r13080 r13442  
    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, lhprdg  ! area etc of new ridges 
    506       REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft  ! area etc of rafted ice 
     505      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, lhprdg, svprdg  ! area etc of new ridges 
     506      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft, svprft  ! area etc of rafted ice 
    507507      ! 
    508508      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges 
     
    579579                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg 
    580580                  lhprdg(ji) = lh_ip_2d(ji,jl1) * afrdg 
     581                  svprdg(ji) = sv_ip_2d(ji,jl1) * afrdg 
    581582                  aprft1     = a_ip_2d(ji,jl1) * afrft 
    582583                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 
    583584                  vprft (ji) = v_ip_2d(ji,jl1) * afrft 
    584585                  lhprft(ji) = lh_ip_2d(ji,jl1) * afrft 
     586                  svprft(ji) = sv_ip_2d(ji,jl1) * afrft 
    585587               ENDIF 
    586588 
     
    613615                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
    614616                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - lhprdg(ji) - lhprft(ji) 
     617                  sv_ip_2d(ji,jl1) = sv_ip_2d(ji,jl1) - svprdg(ji) - svprft(ji) 
    615618               ENDIF 
    616619            ENDIF 
     
    711714                     lh_ip_2d (ji,jl2) = lh_ip_2d(ji,jl2) + (   lhprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
    712715                        &                                   + lhprft (ji) * rn_fpndrft * zswitch(ji)   ) 
     716                     sv_ip_2d (ji,jl2) = sv_ip_2d(ji,jl2) + (   svprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
     717                        &                                   + svprft (ji) * rn_fpndrft * zswitch(ji)   ) 
    713718                  ENDIF 
    714719                   
     
    741746      !---------------- 
    742747      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    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 ) 
     748      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, ,sv_ip_2d, ze_s_2d, ze_i_2d ) 
    744749      ! 
    745750   END SUBROUTINE rdgrft_shift 
     
    854859         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    855860         CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 
     861         CALL tab_3d_2d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip(:,:,:) ) 
    856862         DO jl = 1, jpl 
    857863            DO jk = 1, nlay_s 
     
    868874         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    869875         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) ) 
     876         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_pnd_1d    (1:npti), sfx_pnd    (:,:) ) 
    870877         ! 
    871878         !                 !---------------------! 
     
    881888         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    882889         CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 
     890         CALL tab_2d_3d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip(:,:,:) ) 
    883891         DO jl = 1, jpl 
    884892            DO jk = 1, nlay_s 
     
    895903         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    896904         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) ) 
     905         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_pnd_1d    (1:npti), sfx_pnd    (:,:) ) 
    897906         ! 
    898907      END SELECT 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceistate.F90

    r13080 r13442  
    162162      a_ip_eff (:,:,:) = 0._wp 
    163163      h_ip     (:,:,:) = 0._wp 
     164      s_ip     (:,:,:) = 0._wp 
    164165      ! 
    165166      ! ice velocities 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceitd.F90

    r13080 r13442  
    412412      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    413413      CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 
     414      CALL tab_3d_2d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip ) 
    414415      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
    415416      DO jl = 1, jpl 
     
    488489                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - ztrans 
    489490                  lh_ip_2d(ji,jl2) = lh_ip_2d(ji,jl2) + ztrans 
     491                  !                                               
     492                  ztrans          = sv_ip_2d(ji,jl1) * zworka(ji)     ! Pond salt content 
     493                  sv_ip_2d(ji,jl1) = sv_ip_2d(ji,jl1) - ztrans 
     494                  sv_ip_2d(ji,jl2) = sv_ip_2d(ji,jl2) + ztrans 
    490495               ENDIF 
    491496               ! 
     
    532537      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
    533538      !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
    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 ) 
     539      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, sv_ip_2d, ze_s_2d, ze_i_2d ) 
    535540 
    536541      ! at_i must be <= rn_amax 
     
    561566      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    562567      CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 
     568      CALL tab_2d_3d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip ) 
    563569      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
    564570      DO jl = 1, jpl 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icerst.F90

    r13080 r13442  
    133133      CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip  ) 
    134134      CALL iom_rstput( iter, nitrst, numriw, 'lh_ip', lh_ip ) 
     135      CALL iom_rstput( iter, nitrst, numriw, 'sv_ip', sv_ip  ) 
    135136      ! Snow enthalpy 
    136137      DO jk = 1, nlay_s  
     
    259260            lh_ip(:,:,:) = 0._wp 
    260261         ENDIF 
     262         ! melt pond salinity 
     263         id5 = iom_varid( numrir, 'sv_ip' , ldstop = .FALSE. ) 
     264         IF( id5 > 0 ) THEN 
     265            CALL iom_get( numrir, jpdom_autoglo, 'sv_ip', sv_ip) 
     266         ELSE 
     267            sv_ip(:,:,:) = 0._wp 
     268         ENDIF 
    261269         ! fields needed for Met Office (Jules) coupling 
    262270         IF( ln_cpl ) THEN 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icestp.F90

    r13080 r13442  
    404404      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    405405      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
     406      sfx_pnd(:,:) = 0._wp 
    406407      ! 
    407408      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd.F90

    r13080 r13442  
    217217            ! 
    218218            s_i_new   (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp  ! --- some init --- !  (important to have them here)  
     219            s_i_pndfrz(1:npti) = 0._wp 
    219220            dh_i_sum  (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm  (1:npti) = 0._wp  
    220221            dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 
    221222            dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
     223            dh_s_to_pnd(1:npti) = 0._wp ; dh_i_to_pnd(1:npti) = 0._wp 
    222224            !                                       
    223225                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- ! 
     
    357359         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 
    358360         CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d     (1:npti), lh_ip     (:,:,kl) ) 
     361         CALL tab_2d_1d( npti, nptidx(1:npti), s_ip_1d     (1:npti), s_ip     (:,:,kl) ) 
    359362         ! 
    360363         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            ) 
     
    396399         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub          ) 
    397400         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam          ) 
     401         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd          ) 
    398402         ! 
    399403         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd       ) 
     
    465469         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 
    466470         CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d    (1:npti), lh_ip    (:,:,kl) ) 
     471         CALL tab_1d_2d( npti, nptidx(1:npti), s_ip_1d     (1:npti), s_ip     (:,:,kl) ) 
    467472         ! 
    468473         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
     
    490495         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub        ) 
    491496         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam        ) 
     497         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd        ) 
    492498         ! 
    493499         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd     ) 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_dh.F90

    r11715 r13442  
    2424   USE lib_mpp        ! MPP library 
    2525   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
     26   USE icethd_pnd, ONLY : a_pnd_avail 
    2627    
    2728   IMPLICIT NONE 
     
    9899      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness 
    99100      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah 
     101      REAL(wp)                         ::   zdeltah_to_ocn 
    100102      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting  
    101103 
     
    112114         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile 
    113115      END SELECT 
     116 
     117      ! Initialise fraction of sea ice available for melt ponding  
     118      DO ji = 1, npti  
     119         a_pnd_avail_1d(ji) = a_pnd_avail * at_i_1d(ji)  
     120      END DO  
    114121 
    115122      ! initialize layer thicknesses and enthalpies 
     
    241248               zdeltah  (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change 
    242249               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting 
     250               zdeltah_to_ocn = zdeltah(ji,jk) * (1._wp - a_pnd_avail_1d(ji))           ! Only some of this melt makes it into the ocean.  
     251               dh_s_to_pnd(ji) = dh_s_to_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji)                     ! The rest goes onto melt ponds. 
     252 
    243253               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
    244254                
    245255               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0) 
    246                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
     256               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah_to_ocn * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
    247257                
    248258               ! updates available heat + thickness 
     
    345355                
    346356               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     357               zdeltah_to_ocn = zdeltah(ji,jk) * (1._wp - a_pnd_avail_1d(ji))           ! Only some of this melt makes it into the ocean.  
     358               dh_i_to_pnd(ji) = dh_i_to_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji)                     ! The rest goes onto melt ponds. 
    347359                
    348360               zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 
     
    350362               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt 
    351363                
    352                zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
     364               zfmdt          = - rhoi * zdeltah_to_ocn               ! Recompute mass flux [kg/m2, >0] 
    353365                
    354366               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    355367                
    356                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux >0 
     368               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah_to_ocn * s_i_1d(ji) * r1_rdtice    ! Salt flux >0 
    357369               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok) 
    358370               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0 
    359371               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], > 0   
    360372               !  
    361                wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
     373               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah_to_ocn * r1_rdtice                 ! Mass flux 
    362374                
    363375            END IF 
     
    403415      END DO 
    404416 
     417      ! Melt pond freezing 
     418      !------------------ 
     419      ! Melt pond freezing amounts already calculated in ice_thd_pnd_frz 
     420 
     421      ! update ice enthalpy, thickness and salinity 
     422      DO ji = 1, npti 
     423         s_i_pndfrz(ji) = (h_i_1d(ji) * s_i_1d(ji) + dh_i_from_pnd(ji) + s_ip_1d(ji)) / (h_i_1d(ji) + dh_i_from_pnd(ji)) 
     424         ztmelts = - rTmlt * sz_i_1d(ji,1)                                       ! Freezing point [C] 
     425         zEi           = rcp * ztmelts                                           ! New ice is at freezing point [K] 
     426         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_i_from_pnd(ji) * (-zEi * rhoi)     ! New enthalpy 
     427          
     428                 
     429      END DO 
     430 
    405431 
    406432      ! Ice Basal growth  
     
    540566      ! -------------------------- 
    541567      DO ji = 1, npti 
    542          h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) ) 
     568         h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) + dh_i_from_pnd(ji) ) 
    543569      END DO   
    544570 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_pnd.F90

    r13080 r13442  
    3737   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme 
    3838   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012) 
     39   REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
    3940 
    4041   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp 
     
    132133      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            - 
    133134      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio 
    134       REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
    135135      REAL(wp), PARAMETER ::   pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 
    136136      REAL(wp), PARAMETER ::   pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation 
    137137      ! 
    138       REAL(wp) ::   tot_mlt          ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 
    139       REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding 
     138      REAL(wp) ::   zfr_mlt          ! Total ice and snow surface melt feeding into poinds 
    140139      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change) 
    141140      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
     
    180179            h_ip_1d(ji)      = 0._wp 
    181180            lh_ip_1d(ji)     = 0._wp 
     181            s_ip_1d(ji)      = 0._wp 
    182182            !                                                     !--------------------------------! 
    183183         ELSE                                                     ! Case ice thickness >= rn_himin ! 
     
    190190            ! 
    191191            ! 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  
     192            zdv_mlt = -( dh_i_to_pnd(ji)*rhoi + dh_s_to_pnd(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    199193            ! 
    200194            !--- Pond gowth ---! 
     195            v_ip_old = v_ip_1d(ji) 
    201196            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
    202197            ! 
    203198            !--- Lid shrinking. ---! 
    204199            IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 
     200            ! 
     201            ! --- Pond salinity changes ---! 
     202            If ( v_ip_1d(ji) > epsi20 ) THEN 
     203               s_ip_1d(ji) = ( s_ip_1d(ji) * v_ip_old + s_i_1d(ji) * dh_i_to_pnd(ji)*rhoi* z1_rhow* a_i_1d(ji) ) / v_ip_1d(ji) 
     204            ELSE 
     205               s_ip_1d(ji) = 0._wp 
     206            END IF 
    205207            ! 
    206208            ! melt pond mass flux (<0) 
     
    218220            IF ( ln_pnd_lids ) THEN 
    219221 
    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 
     222               ! Code to use if using surface heat fluxes to form lids 
     223               IF ( ln_pnd_hfx ) THEN 
     224                   
     225                  ! Shrink the pond 
     226                  vol_change = dh_i_from_pnd(ji) * a_i_1d(ji) 
     227                  v_ip_1d(ji) = v_ip_1d(ji) - vol_change 
     228 
     229                  ! Grow the lid 
     230                  lh_ip_1d(ji) = lh_ip_1d(ji) + vol_change / za_ip 
     231 
    240232               ELSE 
    241                   zdh_frz = 0._wp 
    242                END IF 
     233 
     234                 ! Code to use if using melt pond lids 
     235                 IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 
     236                    t_grad = (zTp+rt0) - t_su_1d(ji) 
     237 
     238                    ! The following equation is a rearranged form of: 
     239                    ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 
     240                    ! where: lid_thickness_start = lh_ip_1d(ji) 
     241                    !        lid_thickness_end = lh_ip_end 
     242                    !        omega_dt is a bunch of terms in the equation that do not change 
     243                    ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 
     244                    ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density).  
     245 
     246                    lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
     247                    zdh_frz = lh_ip_end - lh_ip_1d(ji) 
     248 
     249                    ! Pond shrinking 
     250                    v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 
     251 
     252                    ! Lid growing 
     253                    IF ( ln_pnd_lids )  lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 
     254                 ELSE 
     255                    zdh_frz = 0._wp 
     256                 END IF 
     257 
     258              END IF 
    243259 
    244260            ELSE 
     
    249265            ! 
    250266            ! 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  
     267            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 
    252268            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 
    253269            ! 
     
    263279               ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
    264280               IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN 
     281                   v_ip_old = v_ip_1d(ji) 
    265282                   h_ip_1d(ji) = zpnd_aspect * zfr_mlt 
    266283                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
    267284                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
    268285                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     286                   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * r1_rdtice               ! Mass flux to ocean 
     287                   sfx_pnd_1d(ji) = sfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * s_ip_1d(ji) * r1_rdtice ! Salinity flux to ocean 
    269288               ENDIF 
    270289 
    271290               ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    272291               IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 
     292                   v_ip_old = v_ip_1d(ji) 
    273293                   h_ip_1d(ji) = 0.5_wp * h_i_1d(ji) 
    274294                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
    275295                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
    276296                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     297                   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * r1_rdtice               ! Mass flux to ocean 
     298                   sfx_pnd_1d(ji) = sfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * s_ip_1d(ji) * r1_rdtice ! Salinity flux to ocean 
    277299               ENDIF 
    278300 
     
    297319               ! Do the drainage using Darcy's law 
    298320               IF ( perm > 0._wp ) THEN 
     321                   v_ip_old = v_ip_1d(ji) 
    299322                   dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number 
    300323                   dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative 
     
    304327                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
    305328                   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 
     329                   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * r1_rdtice               ! Mass flux to ocean 
     330                   sfx_pnd_1d(ji) = sfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * s_ip_1d(ji) * r1_rdtice ! Salinity flux to ocean 
    317331               ENDIF 
    318332            ENDIF 
     
    322336                a_i_1d(ji) = 0._wp 
    323337                h_s_1d(ji) = 0._wp 
     338            ENDIF 
     339 
     340            ! If there is no melt pond left then set all pond stats to zero 
     341            IF ( v_ip_1d(ji) <= 0._wp ) THEN 
     342               v_ip_1d(ji)      = 0._wp 
     343               a_ip_1d(ji)      = 0._wp 
     344               a_ip_frac_1d(ji) = 0._wp 
     345               h_ip_1d(ji)      = 0._wp 
     346               lh_ip_1d(ji)     = 0._wp 
     347               s_ip_1d(ji)      = 0._wp 
    324348            ENDIF 
    325349 
     
    344368      ! 
    345369   END SUBROUTINE pnd_H12 
     370 
     371   SUBROUTINE ice_thd_pnd_frz 
     372      !!-------------------------------------------------------------------------- 
     373      !!                  ***  ROUTINE ice_thd_pnd_frz   *** 
     374      !! 
     375      !! ** Purpose : calculate the amount of energy released from pond freezing, 
     376      !!              alter the surface heat flux and generate new ice. 
     377      !! 
     378      !!--------------------------------------------------------------------------- 
     379 
     380     REAL(wp), PARAMETER ::   t_max_frz = -5.0_wp    ! Ice surface temperature below which all surface flux goes into freezing 
     381 
     382     INTEGER  ::   ji   ! loop indices 
     383     REAL(wp) ::   vol_pnd                            ! Volume of melt pond per m2 of ice 
     384     REAL(wp) ::   t_based_frac                       ! Fraction contribution to freezing based on temperature 
     385     REAL(wp) ::   p_based_frac                       ! Fraction contribution to freezing based on pond area 
     386     REAL(wp) ::   total_frac                         ! Total fraction of surface heat used in pond freezing 
     387     REAL(wp) ::   frz_whole_pnd                      ! How much ebergy is released by completely freezing the melt pong 
     388 
     389     DO ji = 1, npti 
     390 
     391       top_frz_1d(ji) = 0._wp 
     392       v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)       ! Get the pond volume per m2 of grid box 
     393       a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)   ! Get the pond fraction 
     394       vol_pnd = h_ip_1d(ji) * a_ip_frac_1d(ji)      ! Get the pond volume per m2 of ice 
     395  
     396       ! Only start freezing ponds when the following criteria have been met: 
     397       IF ( qml_ice_1d(ji) < 0.01_wp .AND.         ! There is virtually no top melt 
     398            qcn_ice_1d(ji) < 0._wp  .AND.          ! There is a negative surface heat flux (the surface is cooling down) 
     399            t_su_1d(ji) < (zTp+rt0) .AND.          ! The temperature of the surrounding ice is less then the reference temperature 
     400            vol_pnd > 0._wp ) THEN                 ! There is a pond to freeze 
     401 
     402         ! Calculate the fraction of heat to use based on how cold the surface temperature is 
     403         t_based_frac = (t_su_1d(ji) - (zTp+rt0))/t_max_frz 
     404         IF ( t_based_frac > 1._wp ) t_based_frac = 1._wp 
     405 
     406         ! Calculate the fraction of heat to use based on how much pond there is. 
     407         ! Use a fractional area 0.1 larger as surrounding ice can contribute to pond freezing 
     408         p_based_frac = a_ip_frac_1d(ji) + 0.1_wp 
     409         IF ( p_based_frac > 1._wp ) p_based_frac = 1._wp 
     410 
     411         ! Combine both fractional contributions 
     412         total_frac = p_based_frac * t_based_frac 
     413 
     414         ! What is the total amount of energy released by completely freezing the melt pond 
     415         frz_whole_pnd = vol_pnd * rhow * rLfus 
     416 
     417         ! How much energy is released from freezing 
     418         frz_ip(ji) = -1.0 * qcn_ice_1d(ji) * total_frac 
     419 
     420         ! Make sure we are not freezng more melt pond than there is 
     421         IF ( frz_ip(i) * rdt_ice > frz_whole_pnd ) THEN 
     422           frz_ip(ji) = frz_whole_pnd * r1_rdtice 
     423         END IF 
     424 
     425         ! Alter the surface energy flux to account for some energy from freezing 
     426         qcn_ice_1d(ji) = qcn_ice_1d(ji) + frz_ip(ji) 
     427          
     428         ! How much sea ice growth is caused by freezing melt ponds 
     429         dh_i_from_pnd(ji) = frz_ip(ji) * rdtice * r1_rLfus * r1_rhoi 
     430 
     431       END IF 
     432     END DO 
     433 
     434   END SUBROUTINE ice_thd_pnd_frz 
    346435 
    347436 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_sal.F90

    r11715 r13442  
    8080               zs_i_si = ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice     
    8181               zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth 
     82               zs_i_pnd = (s_i_pndfrz(ji) - s_i_1d(ji) ) * dh_i_from_pnd  (ji) / h_i_1d(ji) ! bottom growth 
    8283               ! Update salinity (nb: salt flux already included in icethd_dh) 
    83                s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si 
     84               s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si + zs_i_pnd 
    8485            ENDIF 
    8586            ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceupdate.F90

    r13080 r13442  
    165165            !------------------------------------------ 
    166166            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   & 
    167                &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) 
     167               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)   & 
     168               &       + sfx_pnd(ji,jj) 
    168169             
    169170            ! Mass of snow and ice per unit area    
     
    209210      IF( iom_use('sfxres'  ) )   CALL iom_put( 'sfxres', sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes 
    210211      IF( iom_use('sfxsub'  ) )   CALL iom_put( 'sfxsub', sfx_sub * 1.e-03 )   ! salt flux from sublimation 
     212      IF( iom_use('sfxsub'  ) )   CALL iom_put( 'sfxpnd', sfx_pnd * 1.e-03 )   ! salt flux from melt ponds 
    211213 
    212214      ! --- mass fluxes [kg/m2/s] --- ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icevar.F90

    r13080 r13442  
    534534               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
    535535               lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj) 
     536               sv_ip (ji,jj,jl) = sv_ip (ji,jj,jl) * zswitch(ji,jj) 
    536537               ! 
    537538            END DO 
     
    556557 
    557558 
    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 ) 
     559   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 
    559560      !!------------------------------------------------------------------- 
    560561      !!                   ***  ROUTINE ice_var_zapneg *** 
     
    572573      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
    573574      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
     575      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_ip     ! melt pond salt content 
    574576      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    575577      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    639641      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    640642      !                                                        but it does not change conservation, so keep it this way is ok 
     643      WHERE( psv_ip (:,:,:) < 0._wp )   psv_ip (:,:,:) = 0._wp 
    641644      WHERE( plh_ip (:,:,:) < 0._wp )   plh_ip (:,:,:) = 0._wp 
    642645      ! 
     
    644647 
    645648 
    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 ) 
     649   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 
    647650      !!------------------------------------------------------------------- 
    648651      !!                   ***  ROUTINE ice_var_roundoff *** 
     
    658661      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
    659662      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
     663      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_ip     ! melt pond salt content 
    660664      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    661665      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    673677         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
    674678         WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 )    plh_ip(1:npti,:)   = 0._wp   ! lh_ip must be >= 0 
     679         WHERE( psv_ip(1:npti,:) < 0._wp .AND. psv_ip(1:npti,:) > -epsi10 )    psv_ip(1:npti,:)   = 0._wp   ! sv_ip must be >= 0 
    675680      ENDIF 
    676681      ! 
     
    792797   !!------------------------------------------------------------------- 
    793798   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    794       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     799      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, psmip,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 
    795800      !!------------------------------------------------------------------- 
    796801      !! ** Purpose :  converting 1-cat ice to 1 ice category 
     
    798803      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    799804      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    800       REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    801       REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     805      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, psmip    ! input  ice/snow temp & sal & ponds 
     806      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip    ! output ice/snow temp & sal & ponds 
    802807      !!------------------------------------------------------------------- 
    803808      ! == thickness and concentration == ! 
     
    813818      pa_ip(:) = patip(:) 
    814819      ph_ip(:) = phtip(:) 
     820      ps_ip (:) = psmip (:) 
    815821       
    816822   END SUBROUTINE ice_var_itd_1c1c 
    817823 
    818824   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    819       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     825      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, psmip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 
    820826      !!------------------------------------------------------------------- 
    821827      !! ** Purpose :  converting N-cat ice to 1 ice category 
     
    823829      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    824830      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    825       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    826       REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     831      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, psmip    ! input  ice/snow temp & sal & ponds 
     832      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip    ! output ice/snow temp & sal & ponds 
    827833      ! 
    828834      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs 
     
    862868      ELSEWHERE                    ;   ph_ip(:) = 0._wp 
    863869      END WHERE 
     870      WHERE( pa_ip(:) /= 0._wp )   ;   ps_ip(:) = SUM( psmip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     871      ELSEWHERE                    ;   ps_ip(:) = 0._wp 
     872      END WHERE 
    864873      ! 
    865874      DEALLOCATE( z1_ai, z1_vi, z1_vs ) 
     
    868877    
    869878   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    870       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     879      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, psmip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 
    871880      !!------------------------------------------------------------------- 
    872881      !! 
     
    890899      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    891900      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    892       REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    893       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     901      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, psmip    ! input  ice/snow temp & sal & ponds 
     902      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip    ! output ice/snow temp & sal & ponds 
    894903      ! 
    895904      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti 
     
    10011010         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
    10021011         END WHERE 
     1012         ps_ip (:,jl) = psmip (:) 
    10031013      END DO 
    10041014      DEALLOCATE( zfra ) 
     
    10071017 
    10081018   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    1009       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     1019      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, psmip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 
    10101020      !!------------------------------------------------------------------- 
    10111021      !! 
     
    10381048      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    10391049      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    1040       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    1041       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     1050      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, psmip    ! input  ice/snow temp & sal & ponds 
     1051      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip    ! output ice/snow temp & sal & ponds 
    10421052      ! 
    10431053      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2 
     
    10681078         pa_ip(:,:) = patip(:,:) 
    10691079         ph_ip(:,:) = phtip(:,:) 
     1080         ps_ip(:,:) = psmip (:,:) 
    10701081         !                              ! ---------------------- ! 
    10711082      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
     
    10731084         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 
    10741085            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
    1075             &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 
    1076             &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  ) 
     1086            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), psmip(:,1), & 
     1087            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ps_ip(:,:)  ) 
    10771088         !                              ! ---------------------- ! 
    10781089      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
     
    10801091         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 
    10811092            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
    1082             &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 
    1083             &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  ) 
     1093            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), psmip(:,:), & 
     1094            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ps_ip(:,1)  ) 
    10841095         !                              ! ----------------------- ! 
    10851096      ELSE                              ! input cat /= output cat ! 
     
    12011212         END DO 
    12021213         ! 
    1203          DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 
     1214         DEALLOCATE( z1_ai, z1_vi, z1_vs ) 
    12041215         ! 
    12051216         ! == ponds == ! 
     
    12231234            END WHERE 
    12241235         END DO 
    1225          DEALLOCATE( zfra ) 
     1236         ztmp(:) = SUM( psmip (:,:) * patip(:,:) * phtip(:,:), dim=2 ) / SUM( phtip(:,:) * patip(:,:), dim=2 ) 
     1237         DO jl = 1, jpl 
     1238            ps_i (:,jl) = ztmp(:) 
     1239         END DO 
     1240         DEALLOCATE( zfra, ztmp ) 
    12261241         ! 
    12271242      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.