Changeset 12382


Ignore:
Timestamp:
2020-02-13T16:36:44+01:00 (12 days ago)
Author:
dancopsey
Message:

Add pond lid code that does most of the science.

Location:
NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icealb.F90

    r10888 r12382  
    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_add_pond_lids_prints/src/ICE/icesbc.F90

    r10888 r12382  
    128128 
    129129      ! --- cloud-sky and overcast-sky ice albedos --- ! 
    130       CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) 
     130      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, lh_ip, zalb_cs, zalb_os ) 
    131131 
    132132      ! albedo depends on cloud fraction because of non-linear spectral effects 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90

    r12379 r12382  
    129129      REAL(wp), PARAMETER ::   zpnd_aspect = 0.174_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 * rdtice / (rLfus * rhow) 
    144154 
    145155      DO ji = 1, npti 
     
    163173            ! 
    164174            ! available meltwater for melt ponding [m, >0] and fraction 
    165             zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
     175            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 
    166176            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    167177            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
     178            actual_mlt = zfr_mlt * zdv_mlt 
    168179            ! 
    169180            !--- Pond gowth ---! 
    170             ! v_ip should never be negative, otherwise code crashes 
    171             v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 
    172             ! 
    173             ! melt pond mass flux (<0) 
    174             IF( zdv_mlt > 0._wp ) THEN 
    175                zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice 
    176                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    177                ! 
    178                ! adjust ice/snow melting flux to balance melt pond flux (>0) 
    179                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 
    180                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    181                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    182             ENDIF 
     181            v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji) 
     182            ! 
     183            !--- Lid shrinking ---! 
     184            lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt 
     185            ! 
    183186            ! 
    184187            !--- Pond contraction (due to refreezing) ---! 
    185             v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 
     188            IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN 
     189               t_grad = freezing_t - t_su_1d(ji) 
     190                
     191               ! The following equation is a rearranged form of: 
     192               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdtice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 
     193               ! where: lid_thickness_start = lh_ip_1d(ji) 
     194               !        lid_thickness_end = lh_ip_end 
     195               !        omega_dt is a bunch of terms in the equation that do not change 
     196               ! 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 
     197               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification. 
     198                              
     199               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
     200               actual_frz = lh_ip_end - lh_ip_1d(ji) 
     201 
     202               ! Pond shrinking 
     203               v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji) 
     204 
     205               ! Lid growing 
     206               lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz 
     207            ELSE 
     208               actual_frz = 0._wp 
     209            END IF 
     210 
     211            ! melt pond mass flux diagnostic (melt only) 
     212            zfac = actual_mlt * a_i_1d(ji * rhow * r1_rdtice 
     213            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
     214            ! 
     215            ! adjust ice/snow melting flux to balance melt pond flux (>0) 
     216            zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 
     217            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
     218            wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
     219            ! 
     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 
    186223            ! 
    187224            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 
     
    190227            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    191228            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 
    192238            ! 
    193239         ENDIF 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90

    r12371 r12382  
    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_add_pond_lids_prints/src/OCE/SBC/sbccpl.F90

    r12374 r12382  
    21222122      INTEGER, INTENT(in) ::   kt 
    21232123      ! 
     2124      REAL(wp), PARAMETER :: pnd_lid_max = 0.015                  !  pond lid thickness above which the ponds disappear from the albedo calculation 
     2125      REAL(wp), PARAMETER :: pnd_lid_min = 0.005                  !  pond lid thickness below which the full pond area is used in the albedo calculation 
     2126                                                                  ! Note: these two variables are mirrored in icealb.F90 (maybe put them in one place...) 
     2127      ! 
    21242128      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    21252129      INTEGER ::   isec, info   ! local integer 
     
    21272131      REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 
    21282132      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4    
     2133      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   lfrac_pnd                 ! The fraction of the meltpond exposed (not inder a frozen lid) 
     2134 
    21292135      !!---------------------------------------------------------------------- 
    21302136      ! 
     
    23412347            SELECT CASE( sn_snd_mpnd%clcat )   
    23422348            CASE( 'yes' )   
    2343                ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) 
     2349 
     2350               ! Calculate how much meltpond is exposed (not under a frozen lid) 
     2351               lfrac_pnd(:,:,1:jpl) = 1.0 
     2352               WHERE( lh_ip(:,:,1:jpl) > pnd_lid_max ) 
     2353                    lfrac_pnd(:,:,1:jpl) = 0.0 
     2354               END WHERE 
     2355               WHERE( lh_ip(:,:,1:jpl) > pnd_lid_min .AND. lh_ip(:,:,1:jpl) <= pnd_lid_max ) 
     2356                    lfrac_pnd(:,:,1:jpl) = ( lh_ip(:,:,1:jpl) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 
     2357               END WHERE 
     2358 
     2359               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) * lfrac_pnd(:,:,1:jpl) 
    23442360               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl) 
    23452361            CASE( 'no' )   
Note: See TracChangeset for help on using the changeset viewer.