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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90 @ 8598

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

second step for making ponds compliant with the new ice
File size: 14.7 KB
Line 
1MODULE icethd_pnd 
2   !!======================================================================
3   !!                     ***  MODULE  icethd_pnd   ***
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   !!   ice_thd_pnd_init      : some initialization and namelist read
15   !!   ice_thd_pnd           : 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   ice_thd_pnd_init    ! routine called by icestp.F90
31   PUBLIC   ice_thd_pnd         ! 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 ice_thd_pnd( kt )
49      !!-------------------------------------------------------------------
50      !!               ***  ROUTINE ice_thd_pnd   ***
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('ice_thd_pnd')
64
65      SELECT CASE ( nice_pnd )
66
67      CASE (np_pndCST)
68         !                             !-------------------------------!
69         CALL pnd_CST                  ! Constant melt ponds           !
70         !                             !-------------------------------!
71      CASE (np_pndH12)
72         !                             !-------------------------------!
73         CALL pnd_H12                  ! Holland et al 2012 melt ponds !
74         !                             !-------------------------------!
75      END SELECT
76
77      IF( nn_timing == 1 )  CALL timing_stop('ice_thd_pnd')
78
79   END SUBROUTINE ice_thd_pnd 
80
81   SUBROUTINE pnd_CST 
82      !!-------------------------------------------------------------------
83      !!                ***  ROUTINE pnd_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(:,:,:) > 0._wp .AND. t_su(:,:,:) >= rt0 ) 
100         a_ip_frac = rn_apnd
101         h_ip      = rn_hpnd   
102         v_ip      = a_ip_frac * a_i * h_ip 
103         a_ip      = a_ip_frac * a_i
104      ELSEWHERE
105         a_ip      = 0._wp
106         h_ip      = 0._wp
107         v_ip      = 0._wp
108         a_ip_frac = 0._wp
109      END WHERE
110     
111   END SUBROUTINE pnd_CST
112
113   SUBROUTINE pnd_H12
114      !!-------------------------------------------------------------------
115      !!                ***  ROUTINE pnd_H12  ***
116      !!
117      !! ** Purpose    : Compute melt pond evolution
118      !!
119      !! ** Method     : Empirical method. A fraction of meltwater is accumulated
120      !!                 in pond volume. It is then released exponentially when
121      !!                 surface is freezing.
122      !!
123      !! ** Tunable parameters : (no real expertise yet, ideas?)
124      !!
125      !! ** Note       : Stolen from CICE for quick test of the melt pond
126      !!                 radiation and freshwater interfaces
127      !!                 Coupling can be radiative AND freshwater
128      !!                 Advection, ridging, rafting are called
129      !!
130      !! ** References : Holland, M. M. et al (J Clim 2012)
131      !!   
132      !!-------------------------------------------------------------------
133     
134      INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting
135      INTEGER, DIMENSION(jpij)         ::   indxj             !
136
137      REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding
138      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zrfrac            ! fraction of available meltwater retained for melt ponding
139
140      REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
141      REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
142      REAL(wp), PARAMETER ::   zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds
143      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
144
145      REAL(wp) ::   zhi               ! dummy ice thickness
146      REAL(wp) ::   zhs               ! dummy snow depth
147      REAL(wp) ::   zTp, z1_Tp        ! reference temperature
148      REAL(wp) ::   zdTs              ! dummy temperature difference
149      REAL(wp) ::   z1_rhofw          ! inverse freshwater density
150      REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
151      REAL(wp) ::   zvpold            ! dummy pond volume
152
153      INTEGER  ::   ji, jj, jl, ij    ! loop indices
154      INTEGER  ::   icells            ! size of dummy array
155      !!-------------------------------------------------------------------
156      z1_rhofw       = 1. / rhofw 
157      z1_zpnd_aspect = 1. / zpnd_aspect
158      zTp            = -2. 
159      z1_Tp          = 1._wp / zTp 
160
161      a_ip_frac(:,:,:) = 0._wp
162      h_ip     (:,:,:) = 0._wp
163
164      !------------------------------------------------------------------
165      ! Available melt water for melt ponding and corresponding fraction
166      !------------------------------------------------------------------
167
168      zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding
169
170      ! NB: zwfx_mlw can be slightly negative for very small values (why?)
171      ! This can in some occasions give negative
172      ! v_ip in the first category, which then gives crazy pond
173      ! fractions and crashes the code as soon as the melt-pond
174      ! radiative coupling is activated
175      ! if we understand and remove why wfx_sum or wfx_snw could be
176      ! negative, then, we can remove the MAX
177      ! NB: I now changed to wfx_snw_sum, this may fix the problem.
178      ! We should check
179
180      zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:) 
181
182      DO jl = 1, jpl   
183
184         !------------------------------------------------------------------------------
185         ! Identify grid cells where ponds should be updated (can probably be improved)
186         !------------------------------------------------------------------------------
187         indxi(:) = 0
188         indxj(:) = 0
189         icells   = 0
190         DO jj = 1, jpj
191            DO ji = 1, jpi
192               IF ( a_i(ji,jj,jl) > epsi10 ) THEN
193                  icells = icells + 1
194                  indxi(icells) = ji
195                  indxj(icells) = jj
196               ENDIF
197            END DO
198         END DO
199
200         DO ij = 1, icells
201
202            ji = indxi(ij)
203            jj = indxj(ij)
204
205            zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
206            zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
207
208            IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin
209
210               a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds
211               v_ip(ji,jj,jl)      = 0._wp
212               a_ip_frac(ji,jj,jl) = 0._wp
213               h_ip(ji,jj,jl)      = 0._wp
214
215               IF( ln_pnd_fwb )   wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)  !--- Give freshwater to the ocean
216
217            ELSE                        !--- Update pond characteristics
218
219               !--- Add retained melt water to melt ponds
220               ! v_ip should never be positive, otherwise code crashes
221               ! MV: as far as I saw, UM5 can create very small negative v_ip values
222               ! hence I added the max, which was not required with Prather (1 yr run)
223               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
224
225               !--- Shrink pond due to refreezing
226               zdTs           = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
227
228               zvpold         = v_ip(ji,jj,jl)
229
230               v_ip(ji,jj,jl) = v_ip(ji,jj,jl) * EXP( zrexp * zdTs * z1_Tp )
231
232               !--- Dump meltwater due to refreezing ( of course this is wrong
233               !--- but this parameterization is too simple )
234               IF ( ln_pnd_fwb ) THEN
235                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
236               ENDIF
237               
238               a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
239               !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
240
241               h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl)
242
243               a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
244
245            ENDIF
246
247         END DO
248
249      END DO
250
251      !--- Remove retained meltwater from surface fluxes
252
253      IF ( ln_pnd_fwb ) THEN
254
255         wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
256         wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )
257
258      ENDIF
259
260   END SUBROUTINE pnd_H12
261
262   SUBROUTINE ice_thd_pnd_init 
263      !!-------------------------------------------------------------------
264      !!                  ***  ROUTINE ice_thd_pnd_init   ***
265      !!
266      !! ** Purpose : Physical constants and parameters linked to melt ponds
267      !!              over sea ice
268      !!
269      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
270      !!               parameter values called at the first timestep (nit000)
271      !!
272      !! ** input   :   Namelist namthd_pnd 
273      !!-------------------------------------------------------------------
274      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read
275      NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
276      !!-------------------------------------------------------------------
277
278      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
279      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
280901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist', lwp )
281
282      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
283      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
284902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist', lwp )
285      IF(lwm) WRITE ( numoni, namthd_pnd )
286     
287      IF(lwp) THEN                        ! control print
288         WRITE(numout,*)
289         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
290         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
291         WRITE(numout,*) '   Namelist namicethd_pnd:'
292         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12
293         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb
294         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST
295         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
296         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
297         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
298      ENDIF
299      !
300      !                             !== set the choice of ice pond scheme ==!
301      ioptio = 0
302                                                            nice_pnd = np_pndNO
303      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
304      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
305      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
306
307      SELECT CASE( nice_pnd )
308      CASE( np_pndNO )         
309         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF
310         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
311      CASE( np_pndCST)
312         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF
313      END SELECT
314      !
315   END SUBROUTINE ice_thd_pnd_init
316   
317#else
318   !!----------------------------------------------------------------------
319   !!   Default option          Empty module          NO ESIM sea-ice model
320   !!----------------------------------------------------------------------
321#endif 
322
323   !!======================================================================
324END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.