Ignore:
Timestamp:
2017-06-06T09:39:43+02:00 (3 years ago)
Author:
vancop
Message:

Melt pond interfaces practically operational

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90

    r8125 r8142  
    66   !! history :       ! Original code by Daniela Flocco and Adrian Turner 
    77   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1) 
    8    !!            2.0  ! 2016    (O. Lecomte, C. Rousset, M. Vancoppenolle) Implementation in NEMO3.6 
     8   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6 
     9   !!                 ! NB: Only lim_mp_cstt and lim_mp_cesm work in this 
     10   !!                       version 
    911   !!---------------------------------------------------------------------- 
    1012#if defined key_lim3 
     
    137139    
    138140 
    139  
    140141   SUBROUTINE lim_mp( kt ) 
    141142      !!------------------------------------------------------------------- 
     
    153154      INTEGER  ::   ji, jj, jl     ! dummy loop indices 
    154155 
    155 !     REAL(wp), POINTER, DIMENSION(:,:) ::   zfsurf  ! surface heat flux(obsolete, should be droped) 
    156       ! 
    157156      !!------------------------------------------------------------------- 
    158  
    159       IF( nn_timing == 1 )  CALL timing_start('limthd') 
    160157 
    161158      IF( nn_timing == 1 )  CALL timing_start('lim_mp') 
     
    180177      END SELECT 
    181178 
    182       ! we should probably not aggregate here since we do it in lim_var_agg 
    183       ! before output, unless we need the total volume and faction else where 
    184  
    185       ! we should also make sure a_ip and v_ip are properly updated at the end 
    186       ! of the routine 
     179      IF( nn_timing == 1 )  CALL timing_stop('lim_mp') 
    187180 
    188181   END SUBROUTINE lim_mp  
     
    237230       !!                 surface is freezing. 
    238231       !! 
    239        !! ** Tunable parameters : (no expertise yet) 
     232       !! ** Tunable parameters : (no real expertise yet, ideas?) 
    240233       !!  
    241234       !! ** Note       : Stolen from CICE for quick test of the melt pond 
     
    272265       !!------------------------------------------------------------------- 
    273266 
    274        CALL wrk_alloc( jpi*jpj, indxi, indxj) 
    275        CALL wrk_alloc( jpi,jpj,     zwfx_mlw ) 
    276        CALL wrk_alloc( jpi,jpj,jpl, zrfrac   ) 
    277  
    278        z1_rhofw       = 1. / rhofw  
    279        z1_zpnd_aspect = 1. / zpnd_aspect 
    280        zTp            = -2.  
    281  
    282        !------------------------------------------------------------------ 
    283        ! Available melt water for melt ponding and corresponding fraction 
    284        !------------------------------------------------------------------ 
    285  
    286        zwfx_mlw(:,:) = wfx_sum(:,:) + wfx_snw(:,:)        ! available meltwater for melt ponding 
    287  
    288        zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)   
    289  
    290        DO jl = 1, jpl    
    291  
    292           ! v_ip(:,:,jl) ! Initialize things 
    293           ! a_ip(:,:,jl) 
    294           ! volpn(:,:) = hpnd(:,:) * apnd(:,:) * aicen(:,:) 
    295  
    296           !------------------------------------------------------------------------------ 
    297           ! Identify grid cells where ponds should be updated (can probably be improved) 
    298           !------------------------------------------------------------------------------ 
    299  
    300           indxi(:) = 0 
    301           indxj(:) = 0 
    302           icells   = 0 
    303  
    304           DO jj = 1, jpj 
    305             DO ji = 1, jpi 
    306                IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
    307                   icells = icells + 1 
    308                   indxi(icells) = ji 
    309                   indxj(icells) = jj 
    310                ENDIF 
    311             END DO                 ! ji 
    312          END DO                    ! jj 
    313  
    314          DO ij = 1, icells 
    315  
    316             ji = indxi(ij) 
    317             jj = indxj(ij) 
    318  
    319             zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    320             zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
    321  
    322             IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin 
    323  
    324                a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds 
    325                v_ip(ji,jj,jl)      = 0._wp 
    326                a_ip_frac(ji,jj,jl) = 0._wp 
    327                h_ip(ji,jj,jl)      = 0._wp 
    328  
    329                IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean 
    330                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)  
    331  
    332  
    333             ELSE                        !--- Update pond characteristics 
    334  
    335                !--- Add retained melt water to melt ponds 
    336                v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice 
    337  
    338                !--- Shrink pond due to refreezing 
    339                zdTs                = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. ) 
    340                 
    341                zvpold              = v_ip(ji,jj,jl) 
    342  
    343                v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) * EXP( zrexp * zdTs / zTp ) 
    344  
    345                !--- Dump meltwater due to refreezing ( of course this is wrong 
    346                !--- but this parameterization is too simple ) 
    347                IF ( ln_pnd_fw ) & 
    348                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice 
    349  
    350                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) ) 
    351  
    352                h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl) 
    353  
    354                a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 
    355  
    356             !----------------------------------------------------------- 
    357             ! Limit pond depth 
    358             !----------------------------------------------------------- 
    359             ! hpondn = min(hpondn, dpthhi*hi) 
    360  
    361             !--- Give freshwater to the ocean ? 
    362  
    363             ENDIF 
    364  
    365           END DO 
    366  
    367        END DO ! jpl 
    368  
    369        !--- Remove retained meltwater from surface fluxes  
    370  
    371        IF ( ln_pnd_fw ) THEN 
    372  
    373            wfx_snw(:,:) = wfx_snw(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
    374  
    375            wfx_sum(:,:) = wfx_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
    376  
    377        ENDIF 
     267        CALL wrk_alloc( jpi*jpj, indxi, indxj) 
     268        CALL wrk_alloc( jpi,jpj,     zwfx_mlw ) 
     269        CALL wrk_alloc( jpi,jpj,jpl, zrfrac   ) 
     270 
     271        z1_rhofw       = 1. / rhofw  
     272        z1_zpnd_aspect = 1. / zpnd_aspect 
     273        zTp            = -2.  
     274 
     275        a_ip_frac(:,:,:) = 0._wp 
     276        h_ip     (:,:,:) = 0._wp 
     277  
     278        !------------------------------------------------------------------ 
     279        ! Available melt water for melt ponding and corresponding fraction 
     280        !------------------------------------------------------------------ 
     281  
     282        zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw(:,:), 0._wp )        ! available meltwater for melt ponding 
     283 
     284                ! NB: zwfx_mlw can be slightly negative for very small values (why?) 
     285                ! This can in some occasions give negative 
     286                ! v_ip in the first category, which then gives crazy pond 
     287                ! fractions and crashes the code as soon as the melt-pond 
     288                ! radiative coupling is activated 
     289                ! if we understand and remove why wfx_sum or wfx_snw could be 
     290                ! negative, then, we can remove the MAX 
     291  
     292        zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)   
     293  
     294        DO jl = 1, jpl    
     295 
     296           !------------------------------------------------------------------------------ 
     297           ! Identify grid cells where ponds should be updated (can probably be improved) 
     298           !------------------------------------------------------------------------------ 
     299  
     300           indxi(:) = 0 
     301           indxj(:) = 0 
     302           icells   = 0 
     303  
     304           DO jj = 1, jpj 
     305             DO ji = 1, jpi 
     306                IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
     307                   icells = icells + 1 
     308                   indxi(icells) = ji 
     309                   indxj(icells) = jj 
     310                ENDIF 
     311             END DO                 ! ji 
     312          END DO                    ! jj 
     313  
     314          DO ij = 1, icells 
     315  
     316             ji = indxi(ij) 
     317             jj = indxj(ij) 
     318  
     319             zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     320             zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
     321  
     322             IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin 
     323  
     324                a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds 
     325                v_ip(ji,jj,jl)      = 0._wp 
     326                a_ip_frac(ji,jj,jl) = 0._wp 
     327                h_ip(ji,jj,jl)      = 0._wp 
     328  
     329                IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean 
     330                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)  
     331  
     332  
     333             ELSE                        !--- Update pond characteristics 
     334  
     335                !--- Add retained melt water to melt ponds 
     336                ! v_ip should never be positive, otherwise code crashes 
     337                ! MV: as far as I saw, UM5 can create very small negative v_ip values 
     338                ! hence I added the max, which was not required with Prather (1 yr run) 
     339                v_ip(ji,jj,jl)      = MAX( v_ip(ji,jj,jl), 0._wp ) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice 
     340  
     341                !--- Shrink pond due to refreezing 
     342                zdTs                = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. ) 
     343                 
     344                zvpold              = v_ip(ji,jj,jl) 
     345  
     346                v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) * EXP( zrexp * zdTs / zTp ) 
     347  
     348                !--- Dump meltwater due to refreezing ( of course this is wrong 
     349                !--- but this parameterization is too simple ) 
     350                IF ( ln_pnd_fw ) & 
     351                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice 
     352  
     353                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) ) 
     354  
     355                h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl) 
     356  
     357                a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 
     358  
     359             !----------------------------------------------------------- 
     360             ! Limit pond depth 
     361             !----------------------------------------------------------- 
     362             ! The original version has pond depth limitation, which I did not 
     363             ! keep here. Maybe we need it later on 
     364             ! 
     365              
     366             ENDIF 
     367  
     368           END DO 
     369  
     370        END DO ! jpl 
     371  
     372        !--- Remove retained meltwater from surface fluxes  
     373  
     374        IF ( ln_pnd_fw ) THEN 
     375  
     376            wfx_snw(:,:) = wfx_snw(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
     377  
     378            wfx_sum(:,:) = wfx_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
     379  
     380        ENDIF 
     381 
     382       CALL wrk_dealloc( jpi*jpj, indxi, indxj) 
     383       CALL wrk_dealloc( jpi,jpj,     zwfx_mlw ) 
     384       CALL wrk_dealloc( jpi,jpj,jpl, zrfrac   ) 
    378385 
    379386   END SUBROUTINE lim_mp_cesm 
Note: See TracChangeset for help on using the changeset viewer.