source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 8370

Last change on this file since 8370 was 8370, checked in by clem, 4 years ago

put limthd_da into the same loop as the other thermodynamics routines

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