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

Changeset 12401


Ignore:
Timestamp:
2020-02-18T16:21:31+01:00 (5 years ago)
Author:
dancopsey
Message:

Add melt pond lid code.

Location:
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src
Files:
18 edited

Legend:

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

    r11715 r12401  
    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   !!-------------|-------------|---------------------------------|-------| 
     
    332333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i) 
    333334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m] 
     335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   lh_ip      !: melt pond lid thickness                  [m] 
    334336 
    335337   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration 
     
    448450 
    449451      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) ) 
     452      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), lh_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 
    451453 
    452454      ii = ii + 1 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice1d.F90

    r11715 r12401  
    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(:) ::   lh_ip_1d      !: Ice pond lid thickness   [m] 
    130131 
    131132   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     
    157158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d 
    158159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d  
     160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lh_ip_2d  
    159161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d  
    160162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_2d 
     
    208210         &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &     
    209211         &      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) ,                   & 
     212         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , lh_ip_1d(jpij) ,  & 
    211213         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   & 
    212214         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
     
    226228      ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  & 
    227229         &      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) ,                      & 
     230         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl),  & 
    229231         &      STAT=ierr(ii) ) 
    230232 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icealb.F90

    r11715 r12401  
    4545CONTAINS 
    4646 
    47    SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os ) 
     47   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, plh_pnd, palb_cs, palb_os ) 
    4848      !!---------------------------------------------------------------------- 
    4949      !!               ***  ROUTINE ice_alb  *** 
     
    9797      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area) 
    9898      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth 
     99      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   plh_pnd      !  melt pond lid thickness 
    99100      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky 
    100101      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_os      !  albedo of ice under overcast sky 
     102      ! 
     103      REAL(wp), PARAMETER :: pnd_lid_max = 0.015                  !  pond lid thickness above which the ponds disappear from the albedo calculation 
     104      REAL(wp), PARAMETER :: pnd_lid_min = 0.005                  !  pond lid thickness below which the full pond area is used in the albedo calculation 
     105                                                                  ! Note: these two variables are mirrored in sbccpl.F90 (maybe put them in one place...) 
    101106      ! 
    102107      INTEGER  ::   ji, jj, jl                ! dummy loop indices 
     
    106111      REAL(wp) ::   zalb_ice, zafrac_ice      ! bare sea ice albedo & relative ice fraction 
    107112      REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction 
     113      REAL(wp) ::   lfrac_pnd                 ! The fraction of the meltpond exposed (not inder a frozen lid) 
    108114      !!--------------------------------------------------------------------- 
    109115      ! 
     
    123129                  zafrac_snw = 0._wp 
    124130                  IF( ld_pnd_alb ) THEN 
    125                      zafrac_pnd = pafrac_pnd(ji,jj,jl) 
     131                     IF ( plh_pnd(ji,jj,jl) > pnd_lid_max ) THEN 
     132                        lfrac_pnd = 0._wp 
     133                     ELSE 
     134                        IF ( plh_pnd(ji,jj,jl) < pnd_lid_min ) THEN 
     135                           lfrac_pnd = 1._wp 
     136                        ELSE 
     137                           lfrac_pnd = ( plh_pnd(ji,jj,jl) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 
     138                        END IF 
     139                     END IF 
     140                     zafrac_pnd = pafrac_pnd(ji,jj,jl) * lfrac_pnd 
    126141                  ELSE 
    127142                     zafrac_pnd = 0._wp 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn.F90

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

    r11715 r12401  
    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_add_pond_lids/src/ICE/icedyn_adv_pra.F90

    r11715 r12401  
    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_add_pond_lids/src/ICE/icedyn_adv_umx.F90

    r11715 r12401  
    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_add_pond_lids/src/ICE/icedyn_rdgrft.F90

    r11715 r12401  
    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_add_pond_lids/src/ICE/iceistate.F90

    r11715 r12401  
    155155      a_ip     (:,:,:) = 0._wp 
    156156      v_ip     (:,:,:) = 0._wp 
     157      lh_ip    (:,:,:) = 0._wp 
    157158      a_ip_frac(:,:,:) = 0._wp 
    158159      h_ip     (:,:,:) = 0._wp 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceitd.F90

    r11715 r12401  
    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_add_pond_lids/src/ICE/icerst.F90

    r11715 r12401  
    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  
     
    245246            CALL iom_get( numrir, jpdom_autoglo, 'a_ip' , a_ip ) 
    246247            CALL iom_get( numrir, jpdom_autoglo, 'v_ip' , v_ip ) 
     248            CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip) 
    247249         ELSE                                     ! start from rest 
    248250            IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it to zero' 
    249251            a_ip(:,:,:) = 0._wp 
    250252            v_ip(:,:,:) = 0._wp 
     253            lh_ip(:,:,:) = 0._wp 
    251254         ENDIF 
    252255         ! fields needed for Met Office (Jules) coupling 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icesbc.F90

    r11715 r12401  
    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_frac, h_ip, lh_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_add_pond_lids/src/ICE/icethd.F90

    r11715 r12401  
    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), lh_ip_1d     (1:npti), lh_ip     (:,:,kl) ) 
    357358         ! 
    358359         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            ) 
     
    461462         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    462463         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
     464         CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d    (1:npti), lh_ip    (:,:,kl) ) 
    463465         ! 
    464466         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd_pnd.F90

    r11715 r12401  
    129129      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio 
    130130      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
     131      REAL(wp), PARAMETER ::   freezing_t  = 273.0_wp ! temperature below which refreezing occurs 
    131132      ! 
    132133      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding 
    133134      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding 
     135      REAL(wp) ::   actual_mlt       ! actual melt water used to make melt ponds (per m2). 
     136                                     ! Need to multiply this by ice area to work on volumes. 
    134137      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    135138      REAL(wp) ::   z1_rhow          ! inverse freshwater density 
    136139      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio 
    137140      REAL(wp) ::   zfac, zdum 
     141      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing 
     142      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing 
     143      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable) 
     144      REAL(wp) ::   actual_frz       ! Amount of melt pond that freezes 
    138145      ! 
    139146      INTEGER  ::   ji   ! loop indices 
     
    142149      z1_zpnd_aspect = 1._wp / zpnd_aspect 
    143150      z1_Tp          = 1._wp / zTp  
     151 
     152      ! Define time-independent field for use in refreezing 
     153      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow) 
    144154 
    145155      DO ji = 1, npti 
     
    151161            a_ip_frac_1d(ji) = 0._wp 
    152162            h_ip_1d(ji)      = 0._wp 
     163            lh_ip_1d(ji)     = 0._wp 
     164 
     165            actual_mlt = 0._wp 
     166            actual_frz = 0._wp 
    153167            !                                                     !--------------------------------! 
    154168         ELSE                                                     ! Case ice thickness >= rn_himin ! 
     
    157171            ! 
    158172            ! 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) 
     173            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 
    160174            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    161175            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
     176            actual_mlt = zfr_mlt * zdv_mlt 
    162177            ! 
    163178            !--- 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 ) 
    166             ! 
     179            v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji) 
     180            ! 
     181            !--- Lid shrinking ---! 
     182            lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt 
     183            ! 
     184            ! 
     185            !--- Pond contraction (due to refreezing) ---! 
     186            IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN 
     187               t_grad = freezing_t - t_su_1d(ji) 
     188                
     189               ! The following equation is a rearranged form of: 
     190               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 
     191               ! where: lid_thickness_start = lh_ip_1d(ji) 
     192               !        lid_thickness_end = lh_ip_end 
     193               !        omega_dt is a bunch of terms in the equation that do not change 
     194               ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same 
     195               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification. 
     196                              
     197               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
     198               actual_frz = lh_ip_end - lh_ip_1d(ji) 
     199 
     200               ! Pond shrinking 
     201               v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji) 
     202 
     203               ! Lid growing 
     204               lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz 
     205            ELSE 
     206               actual_frz = 0._wp 
     207            END IF 
     208 
    167209            ! melt pond mass flux (<0) 
    168210            IF( zdv_mlt > 0._wp ) THEN 
    169                zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice 
     211               zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice 
    170212               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    171213               ! 
     
    176218            ENDIF 
    177219            ! 
    178             !--- 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 ) 
     220            ! Make sure pond volume or lid thickness has not gone negative 
     221            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp  
     222            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 
    180223            ! 
    181224            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 
     
    184227            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    185228            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji) 
     229 
     230            ! If lid thickness is ten times greater than pond thickness then remove pond             
     231            IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     232                a_ip_1d(ji)      = 0._wp 
     233                a_ip_frac_1d(ji) = 0._wp 
     234                h_ip_1d(ji)      = 0._wp 
     235                lh_ip_1d(ji)     = 0._wp 
     236                v_ip_1d(ji)      = 0._wp 
     237            END IF 
    186238            ! 
    187239         ENDIF 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceupdate.F90

    r11715 r12401  
    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_frac, h_ip, lh_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(:,:,:) 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icevar.F90

    r11715 r12401  
    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_add_pond_lids/src/OCE/BDY/bdyice.F90

    r11715 r12401  
    9595            ! exchange 3d arrays 
    9696            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. & 
    97                  &                      , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. & 
     97                 &                      , a_ip, 'T', 1., v_ip, 'T', 1., lh_ip,'T', 1., s_i , 'T', 1., t_su, 'T', 1. & 
    9898                 &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                & 
    9999                 &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      ) 
     
    163163            a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration 
    164164            h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth 
     165            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 
    165166            ! 
    166167            sz_i(ji,jj,:,jl) = s_i(ji,jj,jl) 
     
    170171               a_ip(ji,jj,jl) = 0._wp 
    171172               h_ip(ji,jj,jl) = 0._wp 
     173               lh_ip(ji,jj,jl) = 0._wp 
    172174            ENDIF 
    173175            ! 
     
    231233               a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) 
    232234               h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl) 
     235               lh_ip(ji,jj, jl) = lh_ip(ib,jb, jl) 
    233236               ! 
    234237               sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) 
     
    280283               a_ip(ji,jj,  jl) = 0._wp 
    281284               v_ip(ji,jj,  jl) = 0._wp 
     285               lh_ip(ji,jj, jl) = 0._wp 
    282286               t_su(ji,jj,  jl) = rt0 
    283287               t_s (ji,jj,:,jl) = rt0 
  • NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/SBC/sbccpl.F90

    r11715 r12401  
    20472047      INTEGER, INTENT(in) ::   kt 
    20482048      ! 
     2049      REAL(wp), PARAMETER :: pnd_lid_max = 0.015                  !  pond lid thickness above which the ponds disappear from the albedo calculation 
     2050      REAL(wp), PARAMETER :: pnd_lid_min = 0.005                  !  pond lid thickness below which the full pond area is used in the albedo calculation 
     2051                                                                  ! Note: these two variables are mirrored in icealb.F90 (maybe put them in one place...) 
     2052      ! 
    20492053      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    20502054      INTEGER ::   isec, info   ! local integer 
     
    20522056      REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 
    20532057      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4    
     2058      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   lfrac_pnd                 ! The fraction of the meltpond exposed (not inder a frozen lid) 
     2059 
    20542060      !!---------------------------------------------------------------------- 
    20552061      ! 
     
    22562262            SELECT CASE( sn_snd_mpnd%clcat )   
    22572263            CASE( 'yes' )   
    2258                ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl) 
    2259                ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)   
     2264 
     2265               ! Calculate how much meltpond is exposed (not under a frozen lid) 
     2266               lfrac_pnd(:,:,1:jpl) = 1.0 
     2267               WHERE( lh_ip(:,:,1:jpl) > pnd_lid_max ) 
     2268                    lfrac_pnd(:,:,1:jpl) = 0.0 
     2269               END WHERE 
     2270               WHERE( lh_ip(:,:,1:jpl) > pnd_lid_min .AND. lh_ip(:,:,1:jpl) <= pnd_lid_max ) 
     2271                    lfrac_pnd(:,:,1:jpl) = ( lh_ip(:,:,1:jpl) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 
     2272               END WHERE 
     2273 
     2274               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) * lfrac_pnd(:,:,1:jpl) 
     2275               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl) 
     2276 
    22602277            CASE( 'no' )   
    22612278               ztmp3(:,:,:) = 0.0   
Note: See TracChangeset for help on using the changeset viewer.