source: NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd_dh.F90 @ 11355

Last change on this file since 11355 was 11355, checked in by dancopsey, 2 years ago

Add lots of print statements printing out values of files throughout the ocean_ice timestep.. New output will be in files crash_stats_419.out and crash_stats_68.out.

File size: 36.5 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      INTEGER, DIMENSION(1)  :: imax
108      !!------------------------------------------------------------------
109
110      ! Discriminate between time varying salinity and constant
111      SELECT CASE( nn_icesal )                  ! varying salinity or not
112         CASE( 1 , 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
113         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
114      END SELECT
115
116      ! initialize layer thicknesses and enthalpies
117      h_i_old (1:npti,0:nlay_i+1) = 0._wp
118      eh_i_old(1:npti,0:nlay_i+1) = 0._wp
119      DO jk = 1, nlay_i
120         DO ji = 1, npti
121            h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
122            eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
123         END DO
124      END DO
125      !
126      !                       ! ============================================== !
127      !                       ! Available heat for surface and bottom ablation !
128      !                       ! ============================================== !
129      !
130      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
131         !
132         DO ji = 1, npti
133            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
134         END DO
135         !
136      ELSE
137         !
138         DO ji = 1, npti
139            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
140            qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
141            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
142         END DO
143         !
144      ENDIF
145      !
146      DO ji = 1, npti
147         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 
148         zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rdt_ice )
149      END DO
150
151      ! Ice and snow layer thicknesses
152      !-------------------------------
153      DO jk = 1, nlay_i
154         DO ji = 1, npti
155            zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
156         END DO
157      END DO
158
159      DO jk = 1, nlay_s
160         DO ji = 1, npti
161            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
162         END DO
163      END DO
164
165      !                       ! ============ !
166      !                       !     Snow     !
167      !                       ! ============ !
168      !
169      ! Internal melting
170      ! ----------------
171      ! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
172      DO jk = 1, nlay_s
173         DO ji = 1, npti
174            IF( t_s_1d(ji,jk) > rt0 ) THEN
175               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
176               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux
177               ! updates
178               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk)
179               h_s_1d  (ji)    = h_s_1d(ji) - zh_s(ji,jk)
180               zh_s    (ji,jk) = 0._wp
181               e_s_1d  (ji,jk) = 0._wp
182               t_s_1d  (ji,jk) = rt0
183            END IF
184         END DO
185      END DO
186
187      IF(narea == 419) THEN
188         WRITE(9419,*) 'max wfx_snw_sum_1d icethd_dh after inte5rnal melting = ',MAXVAL(  ABS( wfx_snw_sum_1d(:) )  )
189         imax = MAXLOC( ABS( wfx_snw_sum_1d(:) ) )
190         WRITE(9419,*) 'at location = ',imax
191         WRITE(9419,*) 'zh_s = ',zh_s(imax,:)
192         WRITE(9419,*) 'a_i_1d = ',a_i_1d(imax)
193      ENDIF
194
195      IF(narea == 68) THEN
196         WRITE(968,*) 'max wfx_snw_sum_1d icethd_dh after inte5rnal melting = ',MAXVAL(  ABS( wfx_snw_sum_1d(:) )  )
197         imax = MAXLOC( ABS( wfx_snw_sum_1d(:) ) )
198         WRITE(968,*) 'at location = ',imax
199         WRITE(968,*) 'zh_s = ',zh_s(imax,:)
200         WRITE(968,*) 'a_i_1d = ',a_i_1d(imax)
201      ENDIF
202
203      ! Snow precipitation
204      !-------------------
205      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
206
207      zdeltah(1:npti,:) = 0._wp
208      DO ji = 1, npti
209         IF( sprecip_1d(ji) > 0._wp ) THEN
210            !
211            ! --- precipitation ---
212            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)   ! thickness change
213            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3)
214            !
215            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)
216            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0
217           
218            ! --- melt of falling snow ---
219            rswitch              = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
220            zdeltah       (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change
221            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                 ! bound melting
222            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)
223            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
224           
225            ! updates available heat + precipitations after melting
226            dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)
227            zq_top   (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )     
228            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
229           
230            ! update thickness
231            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
232            !
233         ELSE
234            !
235            zdh_s_pre(ji) = 0._wp
236            zqprec   (ji) = 0._wp
237            !
238         ENDIF
239      END DO
240
241      ! recalculate snow layers
242      DO jk = 1, nlay_s
243         DO ji = 1, npti
244            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
245         END DO
246      END DO
247
248      ! Snow melting
249      ! ------------
250      ! If heat still available (zq_top > 0), then melt more snow
251      zdeltah(1:npti,:) = 0._wp
252      zdh_s_mel(1:npti) = 0._wp
253      DO jk = 1, nlay_s
254         DO ji = 1, npti
255            IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
256               !
257               rswitch          = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) )
258               zdeltah  (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change
259               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting
260               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)
261               
262               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)
263               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)
264               
265               ! updates available heat + thickness
266               dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk)
267               zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
268               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
269               zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
270               !
271            ENDIF
272         END DO
273      END DO
274
275      ! Snow sublimation
276      !-----------------
277      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
278      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
279      zdeltah(1:npti,:) = 0._wp
280      DO ji = 1, npti
281         IF( evap_ice_1d(ji) > 0._wp ) THEN
282            !
283            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice )
284            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)
285            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
286           
287            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)
288               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  &
289               &                 * a_i_1d(ji) * r1_rdtice
290            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
291           
292            ! new snow thickness
293            h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
294            ! update precipitations after sublimation and correct sublimation
295            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
296            zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
297            !
298         ELSE
299            !
300            zdh_s_sub (ji) = 0._wp
301            zevap_rema(ji) = 0._wp
302            !
303         ENDIF
304      END DO
305     
306      ! --- Update snow diags --- !
307      DO ji = 1, npti
308         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
309      END DO
310
311      ! Update temperature, energy
312      !---------------------------
313      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
314      DO jk = 1, nlay_s
315         DO ji = 1,npti
316            rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
317            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            &
318              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  &
319              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) )
320         END DO
321      END DO
322     
323      !                       ! ============ !
324      !                       !     Ice      !
325      !                       ! ============ !
326
327      ! Surface ice melting
328      !--------------------
329      zdeltah(1:npti,:) = 0._wp ! important
330      DO jk = 1, nlay_i
331         DO ji = 1, npti
332            ztmelts = - rTmlt * sz_i_1d(ji,jk)   ! Melting point of layer k [C]
333           
334            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
335
336               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]       
337               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0)
338                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
339               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
340                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
341               zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0
342                         
343               dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting
344               
345               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
346
347               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
348               !                                                                                                  ice enthalpy zEi is "sent" to the ocean
349               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux
350               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
351               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
352
353            ELSE                                        !-- Surface melting
354               
355               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]
356               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
357               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
358               
359               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0]
360               
361               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0]
362               
363               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]
364               
365               zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat
366               
367               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt
368               
369               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
370               
371               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
372               
373               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
374               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok)
375               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0
376               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 
377               !
378               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
379               
380            END IF
381           
382            ! Ice sublimation
383            ! ---------------
384            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi )
385            zdeltah (ji,jk) = zdeltah (ji,jk) + zdum
386            dh_i_sub(ji)    = dh_i_sub(ji)    + zdum
387           
388            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice  ! Salt flux >0
389            !                                                                                          clem: flux is sent to the ocean for simplicity
390            !                                                                                                but salt should remain in the ice except
391            !                                                                                                if all ice is melted. => must be corrected
392            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
393
394            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_rdtice           ! Mass flux > 0
395
396            ! update remaining mass flux
397            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi
398           
399            ! record which layers have disappeared (for bottom melting)
400            !    => icount=0 : no layer has vanished
401            !    => icount=5 : 5 layers have vanished
402            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
403            icount(ji,jk) = NINT( rswitch )
404            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
405                       
406            ! update heat content (J.m-2) and layer thickness
407            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
408            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
409         END DO
410      END DO
411     
412      ! update ice thickness
413      DO ji = 1, npti
414         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) )
415      END DO
416
417      ! remaining "potential" evap is sent to ocean
418      DO ji = 1, npti
419         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)
420      END DO
421
422
423      ! Ice Basal growth
424      !------------------
425      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
426      ! between the inner conductive flux  (qcn_ice_bot), from the open water heat flux
427      ! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
428      ! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice
429
430      ! If salinity varies in time, an iterative procedure is required, because
431      ! the involved quantities are inter-dependent.
432      ! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
433      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
434      ! -> need for an iterative procedure, which converges quickly
435
436      num_iter_max = 1
437      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
438     
439      DO ji = 1, npti
440         IF(  zf_tt(ji) < 0._wp  ) THEN
441            DO iter = 1, num_iter_max   ! iterations
442
443               ! New bottom ice salinity (Cox & Weeks, JGR88 )
444               !--- zswi1  if dh/dt < 2.0e-8
445               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
446               !--- zswi2  if dh/dt > 3.6e-7
447               zgrr     = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_rdtice , epsi10 ) )
448               zswi2    = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
449               zswi12   = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
450               zswi1    = 1. - zswi2 * zswi12
451               zfracs   = MIN( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
452                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
453
454               s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity
455
456               ztmelts       = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C)
457
458               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
459               
460               zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
461                  &            - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts
462
463               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
464
465               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0)
466
467               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
468               
469            END DO
470            ! Contribution to Energy and Salt Fluxes                                   
471            zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0)
472           
473            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
474            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
475           
476            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
477
478            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice                   ! Mass flux, <0
479
480            ! update heat content (J.m-2) and layer thickness
481            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
482            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
483
484         ENDIF
485
486      END DO
487
488      ! Ice Basal melt
489      !---------------
490      zdeltah(1:npti,:) = 0._wp ! important
491      DO jk = nlay_i, 1, -1
492         DO ji = 1, npti
493            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
494
495               ztmelts = - rTmlt * sz_i_1d(ji,jk)  ! Melting point of layer jk (C)
496
497               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
498
499                  zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0)
500                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
501                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted)
502                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
503                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
504
505                  dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk)
506
507                  zfmdt             = - zdeltah(ji,jk) * rhoi      ! Mass flux x time step > 0
508
509                  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
510                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean
511                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux
512                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
513                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
514
515                  ! update heat content (J.m-2) and layer thickness
516                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
517                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
518
519               ELSE                                        !-- Basal melting
520
521                  zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0)
522                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0)
523                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0)
524
525                  zfmdt           = - zq_bot(ji) / zdE                                        ! Mass flux x time step (kg/m2, >0)
526
527                  zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change
528
529                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change
530                 
531                  zq_bot(ji)      = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors
532
533                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt
534
535                  zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0
536
537                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean
538
539                  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 
540                  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 
541
542                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux
543                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok
544                 
545                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
546
547                  ! update heat content (J.m-2) and layer thickness
548                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
549                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
550               ENDIF
551           
552            ENDIF
553         END DO
554      END DO
555
556      ! Update temperature, energy
557      ! --------------------------
558      DO ji = 1, npti
559         h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) )
560      END DO 
561
562      ! If heat still available then melt more snow
563      !-------------------------------------------
564      zdeltah(1:npti,:) = 0._wp ! important
565      DO ji = 1, npti
566         zq_rema (ji)   = zq_top(ji) + zq_bot(ji) 
567         rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
568         rswitch        = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
569         zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
570         zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
571         dh_s_tot(ji)   = dh_s_tot(ji) + zdeltah(ji,1)
572         h_s_1d  (ji)   = h_s_1d  (ji) + zdeltah(ji,1)
573       
574         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2)
575         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)
576         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice       ! Mass flux
577         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
578         !   
579         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
580         qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
581
582         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
583      END DO
584
585      !
586      ! Snow-Ice formation
587      ! ------------------
588      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
589      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
590      z1_rho = 1._wp / ( rhos+rau0-rhoi )
591      DO ji = 1, npti
592         !
593         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho )
594
595         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
596         h_s_1d(ji)    = h_s_1d(ji) - dh_snowice(ji)
597
598         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
599         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0
600         zEw            = rcp * sst_1d(ji)
601         zQm            = zfmdt * zEw 
602         
603         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux
604
605         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux
606
607         ! Case constant salinity in time: virtual salt flux to keep salinity constant
608         IF( nn_icesal /= 2 )  THEN
609            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
610               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice     ! and get  rn_icesal from the ocean
611         ENDIF
612
613         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
614         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice
615         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_rdtice
616
617         ! update heat content (J.m-2) and layer thickness
618         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
619         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
620         
621      END DO
622
623      !
624      ! Update temperature, energy
625      ! --------------------------
626      DO ji = 1, npti
627         rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
628         t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
629      END DO
630
631      DO jk = 1, nlay_s
632         DO ji = 1,npti
633            ! where there is no ice or no snow
634            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) ) ) )
635            ! mass & energy loss to the ocean
636            hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * &
637               &                              ( 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
638            wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * &
639               &                              ( rhos          * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! mass flux
640            ! update energy (mass is updated in the next loop)
641            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
642            ! recalculate t_s_1d from e_s_1d
643            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
644         END DO
645      END DO
646
647      ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
648      WHERE( h_i_1d(1:npti) == 0._wp )   
649         a_i_1d(1:npti) = 0._wp
650         h_s_1d(1:npti) = 0._wp
651      END WHERE
652      !
653   END SUBROUTINE ice_thd_dh
654
655
656   !!--------------------------------------------------------------------------
657   !! INTERFACE ice_thd_snwblow
658   !!
659   !! ** Purpose :   Compute distribution of precip over the ice
660   !!
661   !!                Snow accumulation in one thermodynamic time step
662   !!                snowfall is partitionned between leads and ice.
663   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
664   !!                but because of the winds, more snow falls on leads than on sea ice
665   !!                and a greater fraction (1-at_i)^beta of the total mass of snow
666   !!                (beta < 1) falls in leads.
667   !!                In reality, beta depends on wind speed,
668   !!                and should decrease with increasing wind speed but here, it is
669   !!                considered as a constant. an average value is 0.66
670   !!--------------------------------------------------------------------------
671!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
672   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
673      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
674      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
675      pout = ( 1._wp - ( pin )**rn_blow_s )
676   END SUBROUTINE ice_thd_snwblow_2d
677
678   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
679      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
680      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
681      pout = ( 1._wp - ( pin )**rn_blow_s )
682   END SUBROUTINE ice_thd_snwblow_1d
683
684#else
685   !!----------------------------------------------------------------------
686   !!   Default option                                NO SI3 sea-ice model
687   !!----------------------------------------------------------------------
688#endif
689
690   !!======================================================================
691END MODULE icethd_dh
Note: See TracBrowser for help on using the repository browser.