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 branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90 @ 8879

Last change on this file since 8879 was 8879, checked in by frrh, 6 years ago

Merge in http://fcm3/projects/NEMO.xm/log/branches/UKMO/dev_r8183_ICEMODEL_svn_removed
revisions 8738:8847 inclusive.

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