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 NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90 @ 12382

Last change on this file since 12382 was 12382, checked in by dancopsey, 4 years ago

Add pond lid code that does most of the science.
File size: 15.2 KB
Line 
1MODULE icethd_pnd 
2   !!======================================================================
3   !!                     ***  MODULE  icethd_pnd   ***
4   !!   sea-ice: Melt ponds on top of sea ice
5   !!======================================================================
6   !! history :       !  2012     (O. Lecomte)       Adaptation from Flocco and Turner
7   !!                 !  2017     (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3' :                                     SI3 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 (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
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      !
58      SELECT CASE ( nice_pnd )
59      !
60      CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==!
61         !
62      CASE (np_pndH12)   ;   CALL pnd_H12    !==  Holland et al 2012 melt ponds  ==!
63         !
64      END SELECT
65      !
66   END SUBROUTINE ice_thd_pnd 
67
68
69   SUBROUTINE pnd_CST 
70      !!-------------------------------------------------------------------
71      !!                ***  ROUTINE pnd_CST  ***
72      !!
73      !! ** Purpose :   Compute melt pond evolution
74      !!
75      !! ** Method  :   Melt pond fraction and thickness are prescribed
76      !!                to non-zero values when t_su = 0C
77      !!
78      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
79      !!               
80      !! ** Note   : Coupling with such melt ponds is only radiative
81      !!             Advection, ridging, rafting... are bypassed
82      !!
83      !! ** References : Bush, G.W., and Trump, D.J. (2017)
84      !!-------------------------------------------------------------------
85      INTEGER  ::   ji        ! loop indices
86      !!-------------------------------------------------------------------
87      DO ji = 1, npti
88         !
89         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
90            a_ip_frac_1d(ji) = rn_apnd
91            h_ip_1d(ji)      = rn_hpnd   
92            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
93         ELSE
94            a_ip_frac_1d(ji) = 0._wp
95            h_ip_1d(ji)      = 0._wp   
96            a_ip_1d(ji)      = 0._wp
97         ENDIF
98         !
99      END DO
100      !
101   END SUBROUTINE pnd_CST
102
103
104   SUBROUTINE pnd_H12
105      !!-------------------------------------------------------------------
106      !!                ***  ROUTINE pnd_H12  ***
107      !!
108      !! ** Purpose    : Compute melt pond evolution
109      !!
110      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds
111      !!                 and sent to ocean when surface is freezing
112      !!
113      !!                 pond growth:      Vp = Vp + dVmelt
114      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
115      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
116      !!                    with Tp = -2degC
117      !! 
118      !! ** Tunable parameters : (no real expertise yet, ideas?)
119      !!
120      !! ** Note       : Stolen from CICE for quick test of the melt pond
121      !!                 radiation and freshwater interfaces
122      !!                 Coupling can be radiative AND freshwater
123      !!                 Advection, ridging, rafting are called
124      !!
125      !! ** References : Holland, M. M. et al (J Clim 2012)
126      !!-------------------------------------------------------------------
127      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
128      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            -
129      REAL(wp), PARAMETER ::   zpnd_aspect = 0.174_wp   ! pond aspect ratio
130      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
131      REAL(wp), PARAMETER ::   freezing_t  = 273.0_wp ! temperature below which refreezing occurs
132      !
133      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
134      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
135      REAL(wp) ::   actual_mlt       ! actual melt water used to make melt ponds (per m2).
136                                     ! Need to multiply this by ice area to work on volumes.
137      REAL(wp) ::   z1_Tp            ! inverse reference temperature
138      REAL(wp) ::   z1_rhow          ! inverse freshwater density
139      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
140      REAL(wp) ::   zfac, zdum
141      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing
142      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing
143      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable)
144      REAL(wp) ::   actual_frz       ! Amount of melt pond that freezes
145      !
146      INTEGER  ::   ji   ! loop indices
147      !!-------------------------------------------------------------------
148      z1_rhow        = 1._wp / rhow 
149      z1_zpnd_aspect = 1._wp / zpnd_aspect
150      z1_Tp          = 1._wp / zTp 
151
152      ! Define time-independent field for use in refreezing
153      omega_dt = 2.0_wp * rcnd_i * rdtice / (rLfus * rhow)
154
155      DO ji = 1, npti
156
157         IF (to_print(ji) == 10) THEN
158            write(numout,*)'icethd_pnd_start: h_ip_1d(ji), a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), ' ', a_ip_frac_1d(ji), ' ', a_ip_1d(ji)
159         END IF
160
161         !                                                        !--------------------------------!
162         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
163            !                                                     !--------------------------------!
164            !--- Remove ponds on thin ice
165            a_ip_1d(ji)      = 0._wp
166            a_ip_frac_1d(ji) = 0._wp
167            h_ip_1d(ji)      = 0._wp
168            lh_ip_1d(ji)     = 0._wp
169            !                                                     !--------------------------------!
170         ELSE                                                     ! Case ice thickness >= rn_himin !
171            !                                                     !--------------------------------!
172            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
173            !
174            ! available meltwater for melt ponding [m, >0] and fraction
175            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow
176            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
177            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper
178            actual_mlt = zfr_mlt * zdv_mlt
179            !
180            !--- Pond gowth ---!
181            v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji)
182            !
183            !--- Lid shrinking ---!
184            lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt
185            !
186            !
187            !--- Pond contraction (due to refreezing) ---!
188            IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN
189               t_grad = freezing_t - t_su_1d(ji)
190               
191               ! The following equation is a rearranged form of:
192               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdtice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
193               ! where: lid_thickness_start = lh_ip_1d(ji)
194               !        lid_thickness_end = lh_ip_end
195               !        omega_dt is a bunch of terms in the equation that do not change
196               ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same
197               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification.
198                             
199               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
200               actual_frz = lh_ip_end - lh_ip_1d(ji)
201
202               ! Pond shrinking
203               v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji)
204
205               ! Lid growing
206               lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz
207            ELSE
208               actual_frz = 0._wp
209            END IF
210
211            ! melt pond mass flux diagnostic (melt only)
212            zfac = actual_mlt * a_i_1d(ji * rhow * r1_rdtice
213            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
214            !
215            ! adjust ice/snow melting flux to balance melt pond flux (>0)
216            zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
217            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
218            wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
219            !
220            ! Make sure pond volume or lid thickness has not gone negative
221            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 
222            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp
223            !
224            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
225            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
226            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
227            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
228            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
229
230            ! If lid thickness is ten times greater than pond thickness then remove pond           
231            IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
232                a_ip_1d(ji)      = 0._wp
233                a_ip_frac_1d(ji) = 0._wp
234                h_ip_1d(ji)      = 0._wp
235                lh_ip_1d(ji)     = 0._wp
236                v_ip_1d(ji)      = 0._wp
237            END IF
238            !
239         ENDIF
240
241         IF (to_print(ji) == 10) THEN
242            write(numout,*)'icethd_pnd: h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji)
243            write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdv_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdv_mlt
244            write(numout,*)'icethd_pnd: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi
245         END IF
246
247      END DO
248      !
249   END SUBROUTINE pnd_H12
250
251
252   SUBROUTINE ice_thd_pnd_init 
253      !!-------------------------------------------------------------------
254      !!                  ***  ROUTINE ice_thd_pnd_init   ***
255      !!
256      !! ** Purpose : Physical constants and parameters linked to melt ponds
257      !!              over sea ice
258      !!
259      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
260      !!               parameter values called at the first timestep (nit000)
261      !!
262      !! ** input   :   Namelist namthd_pnd 
263      !!-------------------------------------------------------------------
264      INTEGER  ::   ios, ioptio   ! Local integer
265      !!
266      NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
267      !!-------------------------------------------------------------------
268      !
269      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
270      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
271901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist', lwp )
272      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
273      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
274902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist', lwp )
275      IF(lwm) WRITE ( numoni, namthd_pnd )
276      !
277      IF(lwp) THEN                        ! control print
278         WRITE(numout,*)
279         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
280         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
281         WRITE(numout,*) '   Namelist namicethd_pnd:'
282         WRITE(numout,*) '      Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
283         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
284         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
285         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
286         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
287      ENDIF
288      !
289      !                             !== set the choice of ice pond scheme ==!
290      ioptio = 0
291                                                            nice_pnd = np_pndNO
292      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
293      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
294      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
295      !
296      SELECT CASE( nice_pnd )
297      CASE( np_pndNO )         
298         IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
299      END SELECT
300      !
301   END SUBROUTINE ice_thd_pnd_init
302   
303#else
304   !!----------------------------------------------------------------------
305   !!   Default option          Empty module          NO SI3 sea-ice model
306   !!----------------------------------------------------------------------
307#endif 
308
309   !!======================================================================
310END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.