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 @ 8623

Last change on this file since 8623 was 8623, checked in by clem, 6 years ago

rewrite melt ponds

File size: 14.2 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 NEMO4
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3' :                                     ESIM 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            ! sea-ice: variables
20   USE ice1D          ! sea-ice: thermodynamics variables
21   USE icetab         ! sea-ice: 1D <==> 2D transformation
22   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
26   USE timing         ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   ice_thd_pnd_init    ! routine called by icestp.F90
32   PUBLIC   ice_thd_pnd         ! routine called by icestp.F90
33
34   INTEGER ::              nice_pnd   ! choice of the type of pond scheme
35   !                                               ! associated indices:
36   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme
37   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme
38   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
39
40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
44   !! $Id: icethd_pnd.F90 8420 2017-10-05 13:07:10Z clem $
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE ice_thd_pnd
50      !!-------------------------------------------------------------------
51      !!               ***  ROUTINE ice_thd_pnd   ***
52      !!               
53      !! ** Purpose :   change melt pond fraction
54      !!               
55      !! ** Method  :   brut force
56      !!
57      !! ** Action  : -
58      !!              -
59      !!-------------------------------------------------------------------
60
61      IF( nn_timing == 1 )  CALL timing_start('ice_thd_pnd')
62
63      SELECT CASE ( nice_pnd )
64
65      CASE (np_pndCST)
66         !                             !-------------------------------!
67         CALL pnd_CST                  ! Constant melt ponds           !
68         !                             !-------------------------------!
69      CASE (np_pndH12)
70         !                             !-------------------------------!
71         CALL pnd_H12                  ! Holland et al 2012 melt ponds !
72         !                             !-------------------------------!
73      END SELECT
74
75      IF( nn_timing == 1 )  CALL timing_stop('ice_thd_pnd')
76
77   END SUBROUTINE ice_thd_pnd 
78
79   SUBROUTINE pnd_CST 
80      !!-------------------------------------------------------------------
81      !!                ***  ROUTINE pnd_CST  ***
82      !!
83      !! ** Purpose    : Compute melt pond evolution
84      !!
85      !! ** Method     : Melt pond fraction and thickness are prescribed
86      !!                 to non-zero values when t_su = 0C
87      !!
88      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
89      !!               
90      !! ** Note       : Coupling with such melt ponds is only radiative
91      !!                 Advection, ridging, rafting... are bypassed
92      !!
93      !! ** References : Bush, G.W., and Trump, D.J. (2017)
94      !!   
95      !!-------------------------------------------------------------------
96      INTEGER  ::   ji        ! loop indices
97      !!-------------------------------------------------------------------
98      DO ji = 1, npti
99         
100         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
101            a_ip_frac_1d(ji) = rn_apnd
102            h_ip_1d(ji)      = rn_hpnd   
103            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
104         ELSE
105            a_ip_frac_1d(ji) = 0._wp
106            h_ip_1d(ji)      = 0._wp   
107            a_ip_1d(ji)      = 0._wp
108         ENDIF
109         
110      END DO
111     
112   END SUBROUTINE pnd_CST
113
114   SUBROUTINE pnd_H12
115      !!-------------------------------------------------------------------
116      !!                ***  ROUTINE pnd_H12  ***
117      !!
118      !! ** Purpose    : Compute melt pond evolution
119      !!
120      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds
121      !!                 and sent to ocean when surface is freezing
122      !!
123      !!                 pond growth:      Vp = Vp + dVmelt
124      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
125      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
126      !!                    with Tp = -2degC
127      !!
128      !! ** Tunable parameters : (no real expertise yet, ideas?)
129      !!
130      !! ** Note       : Stolen from CICE for quick test of the melt pond
131      !!                 radiation and freshwater interfaces
132      !!                 Coupling can be radiative AND freshwater
133      !!                 Advection, ridging, rafting are called
134      !!
135      !! ** References : Holland, M. M. et al (J Clim 2012)
136      !!   
137      !!-------------------------------------------------------------------
138      REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
139      REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
140      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
141      REAL(wp), PARAMETER ::   zTp = -2._wp !
142
143      REAL(wp) ::   zfr_mlt            ! fraction of available meltwater retained for melt ponding
144      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
145      REAL(wp) ::   z1_Tp        ! reference temperature
146      REAL(wp) ::   z1_rhofw          ! inverse freshwater density
147      REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
148      REAL(wp) ::   zvpnd            ! dummy pond volume
149      REAL(wp) ::   zfac, zdum
150
151      INTEGER  ::   ji        ! loop indices
152      !!-------------------------------------------------------------------
153      z1_rhofw       = 1. / rhofw 
154      z1_zpnd_aspect = 1. / zpnd_aspect
155      z1_Tp          = 1._wp / zTp 
156
157      DO ji = 1, npti
158         !                                                        !--------------------------------!
159         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
160            !                                                     !--------------------------------!
161            !--- Remove ponds on thin ice and send water into the ocean
162            IF( ln_pnd_fwb ) THEN
163               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + h_ip_1d(ji) * a_ip_1d(ji) * rhofw * r1_rdtice
164            ENDIF
165            a_ip_1d(ji)      = 0._wp
166            a_ip_frac_1d(ji) = 0._wp
167            h_ip_1d(ji)      = 0._wp
168            !                                                     !--------------------------------!
169         ELSE                                                     ! Case ice thickness >= rn_himin !
170            !                                                     !--------------------------------!
171            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
172            zvpnd       = v_ip_1d(ji)
173
174            ! available meltwater for melt ponding [m, >0] and fraction
175            zdv_mlt = -( dh_i_surf(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji)
176            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  !!clem CICE doc
177!!!            zfr_mlt = zrmin + zrmax * a_i_1d(ji)  !!clem Holland paper
178
179            !--- Pond gowth
180            ! v_ip should never be negative, otherwise code crashes
181            ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather)
182            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
183
184            IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN
185               ! melt ponds mass flux (<0)
186               zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice
187               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
188
189               ! adjust ice/snow melting (>0) to take into account ponding
190               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
191               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
192               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
193            ENDIF
194
195            !--- Pond contraction (due to refreezing)
196            v_ip_1d(ji)      = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
197
198            ! --- Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
199            !        h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
200            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
201            !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
202            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
203            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
204
205            !! clem: there is no clever way to conserve mass here...
206           
207!            IF( ln_pnd_fwb ) THEN
208!               ! melt ponds mass flux (>0 when ponds shrink)
209!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw * r1_rdtice
210!               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac!!
211!
212!               ! adjust ice/snow
213!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw / ( a_i_1d(ji)*h_s_1d(ji)*rhosn + a_i_1d(ji)*h_i_1d(ji)*rhoic )!
214!
215!               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + h_s_1d(ji) * zfac * a_i_1d(ji) * rhosn * r1_rdtice
216!               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     + h_i_1d(ji) * zfac * a_i_1d(ji) * rhoic * r1_rdtice
217!               
218!               !!h_s_1d(ji) = h_s_1d(ji) * ( 1._wp + zfac )
219!               !!h_i_1d(ji) = h_i_1d(ji) * ( 1._wp + zfac )!
220!
221!            ENDIF
222
223         ENDIF
224      END DO
225
226   END SUBROUTINE pnd_H12
227
228   SUBROUTINE ice_thd_pnd_init 
229      !!-------------------------------------------------------------------
230      !!                  ***  ROUTINE ice_thd_pnd_init   ***
231      !!
232      !! ** Purpose : Physical constants and parameters linked to melt ponds
233      !!              over sea ice
234      !!
235      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
236      !!               parameter values called at the first timestep (nit000)
237      !!
238      !! ** input   :   Namelist namthd_pnd 
239      !!-------------------------------------------------------------------
240      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read
241      NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
242      !!-------------------------------------------------------------------
243
244      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
245      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
246901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist', lwp )
247
248      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
249      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
250902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist', lwp )
251      IF(lwm) WRITE ( numoni, namthd_pnd )
252     
253      IF(lwp) THEN                        ! control print
254         WRITE(numout,*)
255         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
256         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
257         WRITE(numout,*) '   Namelist namicethd_pnd:'
258         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12
259         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb
260         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST
261         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
262         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
263         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
264      ENDIF
265      !
266      !                             !== set the choice of ice pond scheme ==!
267      ioptio = 0
268                                                            nice_pnd = np_pndNO
269      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
270      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
271      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
272
273      SELECT CASE( nice_pnd )
274      CASE( np_pndNO )         
275         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF
276         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
277      CASE( np_pndCST)
278         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF
279      END SELECT
280      !
281   END SUBROUTINE ice_thd_pnd_init
282   
283#else
284   !!----------------------------------------------------------------------
285   !!   Default option          Empty module          NO ESIM sea-ice model
286   !!----------------------------------------------------------------------
287#endif 
288
289   !!======================================================================
290END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.