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_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90 @ 9012

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

make melt ponds from Holland2012 and associated freshwater flux conservative
File size: 13.0 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   ! reference temperature
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            ! inverse reference temperature
146      REAL(wp) ::   z1_rhofw         ! inverse freshwater density
147      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
148      REAL(wp) ::   zfac, zdum
149
150      INTEGER  ::   ji   ! loop indices
151      !!-------------------------------------------------------------------
152      z1_rhofw       = 1. / rhofw 
153      z1_zpnd_aspect = 1. / zpnd_aspect
154      z1_Tp          = 1._wp / zTp 
155
156      DO ji = 1, npti
157         !                                                        !--------------------------------!
158         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
159            !                                                     !--------------------------------!
160            !--- Remove ponds on thin ice
161            a_ip_1d(ji)      = 0._wp
162            a_ip_frac_1d(ji) = 0._wp
163            h_ip_1d(ji)      = 0._wp
164            !                                                     !--------------------------------!
165         ELSE                                                     ! Case ice thickness >= rn_himin !
166            !                                                     !--------------------------------!
167            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
168
169            ! available meltwater for melt ponding [m, >0] and fraction
170            zdv_mlt = -( dh_i_surf(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji)
171            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
172            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper
173
174            !--- Pond gowth ---!
175            ! v_ip should never be negative, otherwise code crashes
176            ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather)
177            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
178
179            ! melt pond mass flux (<0)
180            IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN
181               zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice
182               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
183
184               ! adjust ice/snow melting flux to balance melt pond flux (>0)
185               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
186               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
187               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
188            ENDIF
189
190            !--- Pond contraction (due to refreezing) ---!
191            v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
192
193            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
194            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
195            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
196            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
197            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
198
199         ENDIF
200      END DO
201
202   END SUBROUTINE pnd_H12
203
204   SUBROUTINE ice_thd_pnd_init 
205      !!-------------------------------------------------------------------
206      !!                  ***  ROUTINE ice_thd_pnd_init   ***
207      !!
208      !! ** Purpose : Physical constants and parameters linked to melt ponds
209      !!              over sea ice
210      !!
211      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
212      !!               parameter values called at the first timestep (nit000)
213      !!
214      !! ** input   :   Namelist namthd_pnd 
215      !!-------------------------------------------------------------------
216      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read
217      NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
218      !!-------------------------------------------------------------------
219
220      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
221      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
222901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist', lwp )
223
224      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
225      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
226902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist', lwp )
227      IF(lwm) WRITE ( numoni, namthd_pnd )
228     
229      IF(lwp) THEN                        ! control print
230         WRITE(numout,*)
231         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
232         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
233         WRITE(numout,*) '   Namelist namicethd_pnd:'
234         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12
235         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb
236         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST
237         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
238         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
239         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
240      ENDIF
241      !
242      !                             !== set the choice of ice pond scheme ==!
243      ioptio = 0
244                                                            nice_pnd = np_pndNO
245      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
246      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
247      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
248
249      SELECT CASE( nice_pnd )
250      CASE( np_pndNO )         
251         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF
252         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
253      CASE( np_pndCST)
254         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF
255      END SELECT
256      !
257   END SUBROUTINE ice_thd_pnd_init
258   
259#else
260   !!----------------------------------------------------------------------
261   !!   Default option          Empty module          NO ESIM sea-ice model
262   !!----------------------------------------------------------------------
263#endif 
264
265   !!======================================================================
266END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.