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
Line 
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)
8   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
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   !!----------------------------------------------------------------------
17   USE phycst           ! physical constants
18   USE dom_oce          ! ocean space and time domain
19   USE ice              ! LIM-3 variables
20   !
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
30   PUBLIC   lim_mp_init    ! routine called by icestp.F90
31   PUBLIC   lim_mp         ! routine called by icestp.F90
32
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
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
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      !!------------------------------------------------------------------------------------
59      INTEGER, INTENT(in) ::   kt    ! number of iteration
60      INTEGER ::   ji, jj, jl     ! dummy loop indices
61      !!-------------------------------------------------------------------
62
63      IF( nn_timing == 1 )  CALL timing_start('lim_mp')
64
65      SELECT CASE ( nice_pnd )
66
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         !                             !------------------------------!
75      END SELECT
76
77      IF( nn_timing == 1 )  CALL timing_stop('lim_mp')
78
79   END SUBROUTINE lim_mp 
80
81   SUBROUTINE lim_mp_CST 
82       !!-------------------------------------------------------------------
83       !!                ***  ROUTINE lim_mp_CST  ***
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       !!-------------------------------------------------------------------
98
99       WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) ) 
100!!       WHERE ( ( a_i > 0._wp ) .AND. ( t_su >= rt0 ) )
101          a_ip_frac = rn_apnd
102          h_ip      = rn_hpnd   
103          v_ip      = a_ip_frac * a_i * h_ip 
104          a_ip      = a_ip_frac * a_i
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
114   END SUBROUTINE lim_mp_CST
115
116   SUBROUTINE lim_mp_H12
117       !!-------------------------------------------------------------------
118       !!                ***  ROUTINE lim_mp_H12  ***
119       !!
120       !! ** Purpose    : Compute melt pond evolution
121       !!
122       !! ** Method     : Empirical method. A fraction of meltwater is accumulated
123       !!                 in pond volume. It is then released exponentially when
124       !!                 surface is freezing.
125       !!
126       !! ** Tunable parameters : (no real expertise yet, ideas?)
127       !!
128       !! ** Note       : Stolen from CICE for quick test of the melt pond
129       !!                 radiation and freshwater interfaces
130       !!                 Coupling can be radiative AND freshwater
131       !!                 Advection, ridging, rafting are called
132       !!
133       !! ** References : Holland, M. M. et al (J Clim 2012)
134       !!   
135       !!-------------------------------------------------------------------
136
137       INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting
138       INTEGER, DIMENSION(jpij)         ::   indxj             !
139
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
142
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
147
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
155
156       INTEGER  ::   ji, jj, jl, ij    ! loop indices
157       INTEGER  ::   icells            ! size of dummy array
158       !!-------------------------------------------------------------------
159        z1_rhofw       = 1. / rhofw 
160        z1_zpnd_aspect = 1. / zpnd_aspect
161        zTp            = -2. 
162
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 
170        zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding
171
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
179                ! NB: I now changed to wfx_snw_sum, this may fix the problem.
180                ! We should check
181 
182        zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:) 
183 
184        DO jl = 1, jpl   
185
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
199             END DO
200          END DO
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 
217                IF ( ln_pnd_fwb ) & !--- Give freshwater to the ocean
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 )
238                IF ( ln_pnd_fwb ) &
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) ) )
242                   !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
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 
256        IF ( ln_pnd_fwb ) THEN
257 
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 ) )
260
261        ENDIF
262
263   END SUBROUTINE lim_mp_H12
264
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      !!-------------------------------------------------------------------
280
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 )
284
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)' )
309
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   
320#else
321   !!----------------------------------------------------------------------
322   !!   Default option          Empty module          NO ESIM sea-ice model
323   !!----------------------------------------------------------------------
324#endif 
325
326   !!======================================================================
327END MODULE limmp 
Note: See TracBrowser for help on using the repository browser.