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.
limmp.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90 @ 8597

Last change on this file since 8597 was 8597, checked in by clem, 7 years ago

first step to make melt ponds compliant with the new code
File size: 14.9 KB
RevLine 
[7293]1MODULE limmp 
2   !!======================================================================
3   !!                     ***  MODULE  limmp   ***
4   !!   Melt ponds
5   !!======================================================================
6   !! history :       ! Original code by Daniela Flocco and Adrian Turner
7   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
[8142]8   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
[7293]9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3' :                                 LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_mp_init      : some initialization and namelist read
15   !!   lim_mp           : main calling routine
16   !!----------------------------------------------------------------------
[7325]17   USE phycst           ! physical constants
18   USE dom_oce          ! ocean space and time domain
[7293]19   USE ice              ! LIM-3 variables
[8378]20   !
[7293]21   USE lbclnk           ! lateral boundary conditions - MPP exchanges
22   USE lib_mpp          ! MPP library
23   USE in_out_manager   ! I/O manager
24   USE lib_fortran      ! glob_sum
25   USE timing           ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
[8321]30   PUBLIC   lim_mp_init    ! routine called by icestp.F90
31   PUBLIC   lim_mp         ! routine called by icestp.F90
[7293]32
[8597]33   INTEGER ::              nice_pnd   ! choice of the type of pond scheme
34   !                                               ! associated indices:
35   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme
36   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme
37   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
38
[7293]39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
43   !! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
[7325]48   SUBROUTINE lim_mp( kt )
49      !!-------------------------------------------------------------------
50      !!               ***  ROUTINE lim_mp   ***
51      !!               
52      !! ** Purpose :   change melt pond fraction
53      !!               
54      !! ** Method  :   brutal force
55      !!
56      !! ** Action  : -
57      !!              -
58      !!------------------------------------------------------------------------------------
[8597]59      INTEGER, INTENT(in) ::   kt    ! number of iteration
60      INTEGER ::   ji, jj, jl     ! dummy loop indices
[7325]61      !!-------------------------------------------------------------------
[7293]62
[7325]63      IF( nn_timing == 1 )  CALL timing_start('lim_mp')
[7293]64
[8597]65      SELECT CASE ( nice_pnd )
[7325]66
[8597]67      CASE (np_pndCST)
68         !                             !------------------------------!
69         CALL lim_mp_CST               ! Constant melt ponds          !
70         !                             !------------------------------!
71      CASE (np_pndH12)
72         !                             !------------------------------!
73         CALL lim_mp_H12               ! Holland et al 2012 melt ponds!
74         !                             !------------------------------!
[8061]75      END SELECT
76
[8142]77      IF( nn_timing == 1 )  CALL timing_stop('lim_mp')
[7325]78
79   END SUBROUTINE lim_mp 
80
[8597]81   SUBROUTINE lim_mp_CST 
[8098]82       !!-------------------------------------------------------------------
[8597]83       !!                ***  ROUTINE lim_mp_CST  ***
[8098]84       !!
85       !! ** Purpose    : Compute melt pond evolution
86       !!
87       !! ** Method     : Melt pond fraction and thickness are prescribed
88       !!                 to non-zero values when t_su = 0C
89       !!
90       !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
91       !!               
92       !! ** Note       : Coupling with such melt ponds is only radiative
93       !!                 Advection, ridging, rafting... are bypassed
94       !!
95       !! ** References : Bush, G.W., and Trump, D.J. (2017)
96       !!   
97       !!-------------------------------------------------------------------
[8061]98
[8098]99       WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) ) 
[8597]100!!       WHERE ( ( a_i > 0._wp ) .AND. ( t_su >= rt0 ) )
[8106]101          a_ip_frac = rn_apnd
[8098]102          h_ip      = rn_hpnd   
[8106]103          v_ip      = a_ip_frac * a_i * h_ip 
104          a_ip      = a_ip_frac * a_i
[8098]105       ELSE WHERE
106          a_ip      = 0._wp
107          h_ip      = 0._wp
108          v_ip      = 0._wp
109          a_ip_frac = 0._wp
110       END WHERE
111
112       wfx_pnd(:,:) = 0._wp
113
[8597]114   END SUBROUTINE lim_mp_CST
[8098]115
[8597]116   SUBROUTINE lim_mp_H12
[8060]117       !!-------------------------------------------------------------------
[8597]118       !!                ***  ROUTINE lim_mp_H12  ***
[8060]119       !!
120       !! ** Purpose    : Compute melt pond evolution
121       !!
[8061]122       !! ** Method     : Empirical method. A fraction of meltwater is accumulated
123       !!                 in pond volume. It is then released exponentially when
[8098]124       !!                 surface is freezing.
[8060]125       !!
[8142]126       !! ** Tunable parameters : (no real expertise yet, ideas?)
[8060]127       !!
128       !! ** Note       : Stolen from CICE for quick test of the melt pond
[8061]129       !!                 radiation and freshwater interfaces
[8098]130       !!                 Coupling can be radiative AND freshwater
131       !!                 Advection, ridging, rafting are called
[8060]132       !!
133       !! ** References : Holland, M. M. et al (J Clim 2012)
134       !!   
135       !!-------------------------------------------------------------------
136
[8597]137       INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting
138       INTEGER, DIMENSION(jpij)         ::   indxj             !
[8060]139
[8597]140       REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding
141       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zrfrac            ! fraction of available meltwater retained for melt ponding
[8060]142
[8597]143       REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
144       REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
145       REAL(wp), PARAMETER ::   zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds
146       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
[8060]147
[8597]148       REAL(wp) ::   zhi               ! dummy ice thickness
149       REAL(wp) ::   zhs               ! dummy snow depth
150       REAL(wp) ::   zTp               ! reference temperature
151       REAL(wp) ::   zdTs              ! dummy temperature difference
152       REAL(wp) ::   z1_rhofw          ! inverse freshwater density
153       REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
154       REAL(wp) ::   zvpold            ! dummy pond volume
[8060]155
[8597]156       INTEGER  ::   ji, jj, jl, ij    ! loop indices
157       INTEGER  ::   icells            ! size of dummy array
[8060]158       !!-------------------------------------------------------------------
[8142]159        z1_rhofw       = 1. / rhofw 
160        z1_zpnd_aspect = 1. / zpnd_aspect
161        zTp            = -2. 
[8060]162
[8142]163        a_ip_frac(:,:,:) = 0._wp
164        h_ip     (:,:,:) = 0._wp
165 
166        !------------------------------------------------------------------
167        ! Available melt water for melt ponding and corresponding fraction
168        !------------------------------------------------------------------
169 
[8179]170        zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding
[8060]171
[8142]172                ! NB: zwfx_mlw can be slightly negative for very small values (why?)
173                ! This can in some occasions give negative
174                ! v_ip in the first category, which then gives crazy pond
175                ! fractions and crashes the code as soon as the melt-pond
176                ! radiative coupling is activated
177                ! if we understand and remove why wfx_sum or wfx_snw could be
178                ! negative, then, we can remove the MAX
[8179]179                ! NB: I now changed to wfx_snw_sum, this may fix the problem.
180                ! We should check
[8142]181 
182        zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:) 
183 
184        DO jl = 1, jpl   
[8060]185
[8142]186           !------------------------------------------------------------------------------
187           ! Identify grid cells where ponds should be updated (can probably be improved)
188           !------------------------------------------------------------------------------
189           indxi(:) = 0
190           indxj(:) = 0
191           icells   = 0
192           DO jj = 1, jpj
193             DO ji = 1, jpi
194                IF ( a_i(ji,jj,jl) > epsi10 ) THEN
195                   icells = icells + 1
196                   indxi(icells) = ji
197                   indxj(icells) = jj
198                ENDIF
[8597]199             END DO
200          END DO
[8142]201 
202          DO ij = 1, icells
203 
204             ji = indxi(ij)
205             jj = indxj(ij)
206 
207             zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
208             zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
209 
210             IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin
211 
212                a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds
213                v_ip(ji,jj,jl)      = 0._wp
214                a_ip_frac(ji,jj,jl) = 0._wp
215                h_ip(ji,jj,jl)      = 0._wp
216 
[8597]217                IF ( ln_pnd_fwb ) & !--- Give freshwater to the ocean
[8142]218                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 
219 
220 
221             ELSE                        !--- Update pond characteristics
222 
223                !--- Add retained melt water to melt ponds
224                ! v_ip should never be positive, otherwise code crashes
225                ! MV: as far as I saw, UM5 can create very small negative v_ip values
226                ! hence I added the max, which was not required with Prather (1 yr run)
227                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
228 
229                !--- Shrink pond due to refreezing
230                zdTs                = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
231               
232                zvpold              = v_ip(ji,jj,jl)
233 
234                v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) * EXP( zrexp * zdTs / zTp )
235 
236                !--- Dump meltwater due to refreezing ( of course this is wrong
237                !--- but this parameterization is too simple )
[8597]238                IF ( ln_pnd_fwb ) &
[8142]239                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
240 
241                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
[8179]242                   !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
[8142]243 
244                h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl)
245 
246                a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
247             
248             ENDIF
249 
250           END DO
251 
252        END DO ! jpl
253 
254        !--- Remove retained meltwater from surface fluxes
255 
[8597]256        IF ( ln_pnd_fwb ) THEN
[8142]257 
[8597]258           wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
259           wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )
[8179]260
[8142]261        ENDIF
[8060]262
[8597]263   END SUBROUTINE lim_mp_H12
[8060]264
[8597]265   SUBROUTINE lim_mp_init 
266      !!-------------------------------------------------------------------
267      !!                  ***  ROUTINE lim_mp_init   ***
268      !!
269      !! ** Purpose : Physical constants and parameters linked to melt ponds
270      !!      over sea ice
271      !!
272      !! ** Method  :  Read the nammp  namelist and check the melt pond 
273      !!       parameter values called at the first timestep (nit000)
274      !!
275      !! ** input   :   Namelist nammp 
276      !!-------------------------------------------------------------------
277      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read
278      NAMELIST/nammp/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
279      !!-------------------------------------------------------------------
[7325]280
[8597]281      REWIND( numnam_ice_ref )              ! Namelist nammp  in reference namelist : Melt Ponds 
282      READ  ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901)
283901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp  in reference namelist', lwp )
[7325]284
[8597]285      REWIND( numnam_ice_cfg )              ! Namelist nammp  in configuration namelist : Melt Ponds
286      READ  ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 )
287902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp )
288      IF(lwm) WRITE ( numoni, nammp )
289     
290      IF(lwp) THEN                        ! control print
291         WRITE(numout,*)
292         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
293         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
294         WRITE(numout,*) '   Namelist namicethd_pnd:'
295         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12
296         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb
297         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST
298         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
299         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
300         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
301      ENDIF
302      !
303      !                             !== set the choice of ice pond scheme ==!
304      ioptio = 0
305                                                            nice_pnd = np_pndNO
306      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
307      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
308      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
[7325]309
[8597]310      SELECT CASE( nice_pnd )
311      CASE( np_pndNO )         
312         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF
313         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
314      CASE( np_pndCST)
315         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF
316      END SELECT
317      !
318   END SUBROUTINE lim_mp_init
319   
[7325]320#else
321   !!----------------------------------------------------------------------
[8597]322   !!   Default option          Empty module          NO ESIM sea-ice model
[7325]323   !!----------------------------------------------------------------------
324#endif 
325
326   !!======================================================================
327END MODULE limmp 
Note: See TracBrowser for help on using the repository browser.