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_dh.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_dh.F90 @ 12369

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

Add print statements
File size: 39.4 KB
Line 
1MODULE icethd_dh
2   !!======================================================================
3   !!                       ***  MODULE icethd_dh ***
4   !!   seaice : thermodynamic growth and melt
5   !!======================================================================
6   !! History :       !  2003-05  (M. Vancoppenolle) Original code in 1D
7   !!                 !  2005-06  (M. Vancoppenolle) 3D version
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_dh        : vertical sea-ice growth and melt
15   !!   ice_thd_snwblow   : distribute snow fall between ice and ocean
16  !!----------------------------------------------------------------------
17   USE dom_oce        ! ocean space and time domain
18   USE phycst         ! physical constants
19   USE ice            ! sea-ice: variables
20   USE ice1D          ! sea-ice: thermodynamics variables
21   USE icethd_sal     ! sea-ice: salinity profiles
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   
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   ice_thd_dh        ! called by ice_thd
31   PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbccpl and here
32
33   INTERFACE ice_thd_snwblow
34      MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d
35   END INTERFACE
36
37   !!----------------------------------------------------------------------
38   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
39   !! $Id$
40   !! Software governed by the CeCILL license (see ./LICENSE)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE ice_thd_dh
45      !!------------------------------------------------------------------
46      !!                ***  ROUTINE ice_thd_dh  ***
47      !!
48      !! ** Purpose :   compute ice and snow thickness changes due to growth/melting
49      !!
50      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
51      !!                Bottom accretion/ablation arises from flux budget
52      !!                Snow thickness can increase by precipitation and decrease by sublimation
53      !!                If snow load excesses Archmiede limit, snow-ice is formed by
54      !!                the flooding of sea-water in the snow
55      !!
56      !!                - Compute available flux of heat for surface ablation
57      !!                - Compute snow and sea ice enthalpies
58      !!                - Surface ablation and sublimation
59      !!                - Bottom accretion/ablation
60      !!                - Snow ice formation
61      !!
62      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
63      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
64      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
65      !!              Vancoppenolle et al.,2009, Ocean Modelling
66      !!------------------------------------------------------------------
67      INTEGER  ::   ji, jk       ! dummy loop indices
68      INTEGER  ::   iter         ! local integer
69
70      REAL(wp) ::   ztmelts      ! local scalar
71      REAL(wp) ::   zdum       
72      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
73      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
74      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
75      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
76      REAL(wp) ::   zgrr         ! bottom growth rate
77      REAL(wp) ::   zt_i_new     ! bottom formation temperature
78      REAL(wp) ::   z1_rho       ! 1/(rhos+rau0-rhoi)
79
80      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
81      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
82      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
83      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
84      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
85
86      REAL(wp), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3)
87      REAL(wp), DIMENSION(jpij) ::   zq_top      ! heat for surface ablation                   (J.m-2)
88      REAL(wp), DIMENSION(jpij) ::   zq_bot      ! heat for bottom ablation                    (J.m-2)
89      REAL(wp), DIMENSION(jpij) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
90      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
91      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
92
93      REAL(wp), DIMENSION(jpij) ::   zdh_s_mel   ! snow melt
94      REAL(wp), DIMENSION(jpij) ::   zdh_s_pre   ! snow precipitation
95      REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation
96
97      REAL(wp), DIMENSION(jpij,nlay_s) ::   zh_s      ! snw layer thickness
98      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
99      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
100      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting
101
102      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
103
104      REAL(wp) ::   zswitch_sal
105
106      INTEGER  ::   num_iter_max      ! Heat conservation
107      !!------------------------------------------------------------------
108
109      ! Discriminate between time varying salinity and constant
110      SELECT CASE( nn_icesal )                  ! varying salinity or not
111         CASE( 1 , 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
112         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
113      END SELECT
114
115      ! initialize layer thicknesses and enthalpies
116      h_i_old (1:npti,0:nlay_i+1) = 0._wp
117      eh_i_old(1:npti,0:nlay_i+1) = 0._wp
118      DO jk = 1, nlay_i
119         DO ji = 1, npti
120            h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
121            eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
122         END DO
123      END DO
124      !
125      !                       ! ============================================== !
126      !                       ! Available heat for surface and bottom ablation !
127      !                       ! ============================================== !
128      !
129      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
130         !
131         DO ji = 1, npti
132            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
133
134            IF ( to_print(ji) == 10 ) THEN
135               write(numout,*)'loading flux: zq_top(ji), qml_ice_1d(ji) = ',zq_top(ji), ' ', qml_ice_1d(ji)
136            ENDIF
137         END DO
138         !
139      ELSE
140         !
141         DO ji = 1, npti
142            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
143            qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
144            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
145
146            IF ( to_print(ji) == 10 ) THEN
147               write(numout,*)'icethd_dh: qns_ice_1d(ji), qsr_ice_1d(ji), qtr_ice_top_1d(ji), qcn_ice_top_1d(ji) = ',qns_ice_1d(ji), ' ', qsr_ice_1d(ji), ' ', qtr_ice_top_1d(ji), ' ', qcn_ice_top_1d(ji)
148               write(numout,*)'icethd_dh: zdum, t_su_1d(ji), qml_ice_1d(ji), zq_top(ji) = ',zdum, ' ', t_su_1d(ji), qml_ice_1d(ji), zq_top(ji)
149            ENDIF
150
151         END DO
152         !
153      ENDIF
154      !
155      DO ji = 1, npti
156         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 
157         zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rdt_ice )
158         IF ( to_print(ji) == 1 ) THEN
159            write(numout,*)'icethd_dh: qcn_ice_bot_1d(ji), qsb_ice_bot_1d(ji) = ',qcn_ice_bot_1d(ji), qsb_ice_bot_1d(ji)
160            write(numout,*)'icethd_dh: fhld_1d(ji), zq_bot(ji) = ',fhld_1d(ji), ', ', zq_bot(ji)
161         ENDIF
162      END DO
163
164      ! Ice and snow layer thicknesses
165      !-------------------------------
166      DO jk = 1, nlay_i
167         DO ji = 1, npti
168            zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
169         END DO
170      END DO
171
172      DO jk = 1, nlay_s
173         DO ji = 1, npti
174            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
175         END DO
176      END DO
177
178      !                       ! ============ !
179      !                       !     Snow     !
180      !                       ! ============ !
181      !
182      ! Internal melting
183      ! ----------------
184      ! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
185      DO jk = 1, nlay_s
186         DO ji = 1, npti
187            IF( t_s_1d(ji,jk) > rt0 ) THEN
188               hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! heat flux to the ocean [W.m-2], < 0
189               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux
190               ! updates
191               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk)
192               h_s_1d  (ji)    = h_s_1d(ji) - zh_s(ji,jk)
193               zh_s    (ji,jk) = 0._wp
194               e_s_1d  (ji,jk) = 0._wp
195               t_s_1d  (ji,jk) = rt0
196            END IF
197
198           IF (to_print(ji) == 10) THEN
199                write(numout,*)'icethd_dh 8: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
200            END IF
201         END DO
202      END DO
203
204      DO ji = 1, npti
205         IF (to_print(ji) == 10) THEN
206            write(numout,*)'icethd_dh internal_melting: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
207         END IF
208      END DO       
209
210      ! Snow precipitation
211      !-------------------
212      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
213
214      zdeltah(1:npti,:) = 0._wp
215      DO ji = 1, npti
216         IF( sprecip_1d(ji) > 0._wp ) THEN
217            !
218            ! --- precipitation ---
219            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)   ! thickness change
220            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3)
221            !
222            hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat flux from snow precip (>0, W.m-2)
223            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0
224           
225            ! --- melt of falling snow ---
226            rswitch              = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
227            zdeltah       (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change
228            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                 ! bound melting
229            hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat used to melt snow (W.m-2, >0)
230            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip), >0
231           
232            ! updates available heat + precipitations after melting
233            dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)
234            zq_top   (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )     
235            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
236
237            IF ( to_print(ji) == 10 ) THEN
238               write(numout,*)'after precip: zq_top(ji), zdeltah(ji,1), zqprec(ji) = ',zq_top(ji), zdeltah(ji,1), zqprec(ji)
239            ENDIF
240           
241            ! update thickness
242            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
243            !
244         ELSE
245            !
246            zdh_s_pre(ji) = 0._wp
247            zqprec   (ji) = 0._wp
248            !
249         ENDIF
250
251         IF (to_print(ji) == 10) THEN
252            write(numout,*)'icethd_dh snow_precip: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
253         END IF
254
255      END DO
256
257      ! recalculate snow layers
258      DO jk = 1, nlay_s
259         DO ji = 1, npti
260            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
261         END DO
262      END DO
263
264      ! Snow melting
265      ! ------------
266      ! If heat still available (zq_top > 0), then melt more snow
267      zdeltah(1:npti,:) = 0._wp
268      zdh_s_mel(1:npti) = 0._wp
269      DO jk = 1, nlay_s
270         DO ji = 1, npti
271            IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
272               !
273               rswitch          = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) )
274               zdeltah  (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change
275               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting
276               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)
277               
278               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0)
279               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip)
280               
281               ! updates available heat + thickness
282               dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk)
283               zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
284               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
285               zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
286               !
287            ENDIF
288
289            IF ( to_print(ji) == 10 ) THEN
290               write(numout,*)'after snow melt: zq_top(ji), zdeltah(ji,jk), e_s_1d(ji,jk), test = ',zq_top(ji), zdeltah(ji,jk), e_s_1d(ji,jk), ( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp )
291            ENDIF
292
293            IF (to_print(ji) == 10) THEN
294               write(numout,*)'icethd_dh loop: wfx_snw_sum_1d(ji), jk = ',wfx_snw_sum_1d(ji), ' ', jk
295            END IF
296         END DO
297      END DO
298
299      DO ji = 1, npti
300         IF (to_print(ji) == 10) THEN
301            write(numout,*)'icethd_dh snow_ground: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
302         END IF
303      END DO
304
305      ! Snow sublimation
306      !-----------------
307      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
308      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
309      zdeltah(1:npti,:) = 0._wp
310      DO ji = 1, npti
311         IF( evap_ice_1d(ji) > 0._wp ) THEN
312            !
313            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice )
314            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos   ! remaining evap in kg.m-2 (used for ice melting later on)
315            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
316           
317            hfx_sub_1d    (ji) = hfx_sub_1d(ji) + &   ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow)
318               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  &
319               &                 * a_i_1d(ji) * r1_rdtice
320            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice   ! Mass flux by sublimation
321           
322            ! new snow thickness
323            h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
324            ! update precipitations after sublimation and correct sublimation
325            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
326            zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
327            !
328         ELSE
329            !
330            zdh_s_sub (ji) = 0._wp
331            zevap_rema(ji) = 0._wp
332            !
333         ENDIF
334      END DO
335     
336      ! --- Update snow diags --- !
337      DO ji = 1, npti
338         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
339      END DO
340
341      ! Update temperature, energy
342      !---------------------------
343      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
344      DO jk = 1, nlay_s
345         DO ji = 1,npti
346            rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
347            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            &
348              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  &
349              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) )
350
351            IF (to_print(ji) == 10) THEN
352                write(numout,*)'icethd_dh 9: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
353                write(numout,*)'icethd_dh 9: zdh_s_pre(ji), zqprec(ji), h_s_1d(ji) = ',zdh_s_pre(ji), ' ', zqprec(ji), ' ', h_s_1d(ji)
354            END IF
355         END DO
356      END DO
357     
358      !                       ! ============ !
359      !                       !     Ice      !
360      !                       ! ============ !
361
362      ! Surface ice melting
363      !--------------------
364      zdeltah(1:npti,:) = 0._wp ! important
365      DO jk = 1, nlay_i
366         DO ji = 1, npti
367            ztmelts = - rTmlt * sz_i_1d(ji,jk)   ! Melting point of layer k [C]
368           
369            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
370
371               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]       
372               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0)
373                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
374               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
375                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
376               zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0
377                         
378               dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting
379               
380               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
381
382               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0
383               !                                                                                                  ice enthalpy zEi is "sent" to the ocean
384               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux
385               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
386               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
387
388            ELSE                                        !-- Surface melting
389               
390               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]
391               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
392               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
393               
394               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0]
395               
396               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0]
397               
398               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
399               
400               zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat
401               
402               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt
403               
404               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
405               
406               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
407               
408               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux >0
409               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok)
410               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0
411               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], > 0 
412               !
413               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
414               
415            END IF
416
417            IF ( to_print(ji) == 10 ) THEN
418               write(numout,*)'after top melt: zq_top(ji), test = ',zq_top(ji), ( t_i_1d(ji,jk) >= (ztmelts+rt0) ) 
419               write(numout,*)'after top melt: wfx_sum_1d(ji), zdeltah(ji,jk), zfmdt, zdE = ', wfx_sum_1d(ji), zdeltah(ji,jk), zfmdt, zdE
420            ENDIF
421           
422            ! Ice sublimation
423            ! ---------------
424            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi )
425            zdeltah (ji,jk) = zdeltah (ji,jk) + zdum
426            dh_i_sub(ji)    = dh_i_sub(ji)    + zdum
427           
428            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice  ! Salt flux >0
429            !                                                                                          clem: flux is sent to the ocean for simplicity
430            !                                                                                                but salt should remain in the ice except
431            !                                                                                                if all ice is melted. => must be corrected
432            hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice      ! Heat flux [W.m-2], < 0
433
434            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_rdtice           ! Mass flux > 0
435
436            ! update remaining mass flux
437            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi
438           
439            ! record which layers have disappeared (for bottom melting)
440            !    => icount=0 : no layer has vanished
441            !    => icount=5 : 5 layers have vanished
442            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
443            icount(ji,jk) = NINT( rswitch )
444            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
445                       
446            ! update heat content (J.m-2) and layer thickness
447            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
448            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
449         END DO
450      END DO
451     
452      ! update ice thickness
453      DO ji = 1, npti
454         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) )
455      END DO
456
457      ! remaining "potential" evap is sent to ocean
458      DO ji = 1, npti
459         wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
460      END DO
461
462
463      ! Ice Basal growth
464      !------------------
465      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
466      ! between the inner conductive flux  (qcn_ice_bot), from the open water heat flux
467      ! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
468      ! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice
469
470      ! If salinity varies in time, an iterative procedure is required, because
471      ! the involved quantities are inter-dependent.
472      ! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
473      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
474      ! -> need for an iterative procedure, which converges quickly
475
476      num_iter_max = 1
477      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
478     
479      DO ji = 1, npti
480         IF(  zf_tt(ji) < 0._wp  ) THEN
481            DO iter = 1, num_iter_max   ! iterations
482
483               ! New bottom ice salinity (Cox & Weeks, JGR88 )
484               !--- zswi1  if dh/dt < 2.0e-8
485               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
486               !--- zswi2  if dh/dt > 3.6e-7
487               zgrr     = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_rdtice , epsi10 ) )
488               zswi2    = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
489               zswi12   = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
490               zswi1    = 1. - zswi2 * zswi12
491               zfracs   = MIN( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
492                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
493
494               s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity
495
496               ztmelts       = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C)
497
498               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
499               
500               zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
501                  &            - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts
502
503               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
504
505               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0)
506
507               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
508               
509            END DO
510            ! Contribution to Energy and Salt Fluxes                                   
511            zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0)
512           
513            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], >0
514            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], <0
515           
516            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice     ! Salt flux, <0
517
518            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice                   ! Mass flux, <0
519
520            ! update heat content (J.m-2) and layer thickness
521            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
522            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
523
524         ENDIF
525
526      END DO
527
528      ! Ice Basal melt
529      !---------------
530      zdeltah(1:npti,:) = 0._wp ! important
531      DO jk = nlay_i, 1, -1
532         DO ji = 1, npti
533            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
534
535               ztmelts = - rTmlt * sz_i_1d(ji,jk)  ! Melting point of layer jk (C)
536
537               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
538
539                  zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0)
540                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
541                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted)
542                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
543                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
544
545                  dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk)
546
547                  zfmdt             = - zdeltah(ji,jk) * rhoi      ! Mass flux x time step > 0
548
549                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0
550                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean
551                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux
552                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
553                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
554
555                  ! update heat content (J.m-2) and layer thickness
556                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
557                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
558
559               ELSE                                        !-- Basal melting
560
561                  zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0)
562                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0)
563                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0)
564
565                  zfmdt           = - zq_bot(ji) / zdE                                        ! Mass flux x time step (kg/m2, >0)
566
567                  zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change
568
569                  IF ( to_print(ji) == 1 ) THEN
570                    write(numout,*)'icethd_dh: zfmdt, zq_bot(ji), zdE, zdeltah(ji,jk) = ',zfmdt,', ', zq_bot(ji), zdE, zdeltah(ji,jk)
571                  ENDIF
572
573                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change
574                 
575                  zq_bot(ji)      = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors
576
577                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt
578
579                  zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0
580
581                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean
582
583                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
584                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat used in this process [W.m-2], >0 
585
586                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux
587                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok
588                 
589                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
590
591                  IF (to_print(ji) == 1) THEN
592                     write(numout,*)'icethd_dh: wfx_bom_1d(ji), in meters, at_i_1d(ji) = ',wfx_bom_1d(ji),wfx_bom_1d(ji)*rdt_ice/(rhoi*at_i_1d(ji)), at_i_1d(ji)
593                  ENDIF
594
595                  ! update heat content (J.m-2) and layer thickness
596                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
597                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
598               ENDIF
599           
600            ENDIF
601         END DO
602      END DO
603
604     
605
606      ! Update temperature, energy
607      ! --------------------------
608      DO ji = 1, npti
609         h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) )
610      END DO 
611
612      ! If heat still available then melt more snow
613      !-------------------------------------------
614      zdeltah(1:npti,:) = 0._wp ! important
615      DO ji = 1, npti
616         zq_rema (ji)   = zq_top(ji) + zq_bot(ji) 
617         rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
618         rswitch        = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
619         zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
620         zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
621         dh_s_tot(ji)   = dh_s_tot(ji) + zdeltah(ji,1)
622         h_s_1d  (ji)   = h_s_1d  (ji) + zdeltah(ji,1)
623       
624         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2)
625         hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice   ! Heat used to melt snow, W.m-2 (>0)
626         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice       ! Mass flux
627         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
628         !   
629         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
630         qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
631
632         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
633
634         IF (to_print(ji) == 10) THEN
635            write(numout,*)'icethd_dh snow_extra: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
636         END IF
637
638      END DO
639
640      !
641      ! Snow-Ice formation
642      ! ------------------
643      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
644      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
645      z1_rho = 1._wp / ( rhos+rau0-rhoi )
646      DO ji = 1, npti
647         !
648         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho )
649
650         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
651         h_s_1d(ji)    = h_s_1d(ji) - dh_snowice(ji)
652
653         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
654         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0
655         zEw            = rcp * sst_1d(ji)
656         zQm            = zfmdt * zEw 
657         
658         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux
659
660         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux
661
662         ! Case constant salinity in time: virtual salt flux to keep salinity constant
663         IF( nn_icesal /= 2 )  THEN
664            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean
665               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice     ! and get  rn_icesal from the ocean
666         ENDIF
667
668         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
669         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice
670         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_rdtice
671
672         ! update heat content (J.m-2) and layer thickness
673         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
674         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
675         
676      END DO
677
678      !
679      ! Update temperature, energy
680      ! --------------------------
681      DO ji = 1, npti
682         rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
683         t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
684
685         IF (to_print(ji) == 10) THEN
686            write(numout,*)'icethd_dh: t_su_1d(ji), t_i_1d(ji,1) = ',t_su_1d(ji), ' ', t_i_1d(ji,1)
687         END IF
688      END DO
689
690      DO jk = 1, nlay_s
691         DO ji = 1,npti
692            ! where there is no ice or no snow
693            rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) )
694            ! mass & energy loss to the ocean
695            hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * &
696               &                              ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! heat flux to the ocean [W.m-2], < 0
697            wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * &
698               &                              ( rhos          * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! mass flux
699            ! update energy (mass is updated in the next loop)
700            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
701            ! recalculate t_s_1d from e_s_1d
702            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
703
704            IF (to_print(ji) == 10) THEN
705                write(numout,*)'icethd_dh 10: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
706            END IF
707         END DO
708      END DO
709
710      ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
711      WHERE( h_i_1d(1:npti) == 0._wp )   
712         a_i_1d(1:npti) = 0._wp
713         h_s_1d(1:npti) = 0._wp
714      END WHERE
715      !
716   END SUBROUTINE ice_thd_dh
717
718
719   !!--------------------------------------------------------------------------
720   !! INTERFACE ice_thd_snwblow
721   !!
722   !! ** Purpose :   Compute distribution of precip over the ice
723   !!
724   !!                Snow accumulation in one thermodynamic time step
725   !!                snowfall is partitionned between leads and ice.
726   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
727   !!                but because of the winds, more snow falls on leads than on sea ice
728   !!                and a greater fraction (1-at_i)^beta of the total mass of snow
729   !!                (beta < 1) falls in leads.
730   !!                In reality, beta depends on wind speed,
731   !!                and should decrease with increasing wind speed but here, it is
732   !!                considered as a constant. an average value is 0.66
733   !!--------------------------------------------------------------------------
734!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
735   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
736      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
737      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
738      pout = ( 1._wp - ( pin )**rn_blow_s )
739   END SUBROUTINE ice_thd_snwblow_2d
740
741   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
742      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
743      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
744      pout = ( 1._wp - ( pin )**rn_blow_s )
745   END SUBROUTINE ice_thd_snwblow_1d
746
747#else
748   !!----------------------------------------------------------------------
749   !!   Default option                                NO SI3 sea-ice model
750   !!----------------------------------------------------------------------
751#endif
752
753   !!======================================================================
754END MODULE icethd_dh
Note: See TracBrowser for help on using the repository browser.