source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 7256

Last change on this file since 7256 was 7256, checked in by cbricaud, 4 years ago

phaze NEMO routines in CRS branch with nemo_v3_6_STABLE branch at rev 7213 (09-09-2016) (merge -r 5519:7213 )

  • Property svn:keywords set to Id
File size: 37.9 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 sbc_oce        ! Surface boundary condition: ocean fields
21   USE ice            ! LIM variables
22   USE thd_ice        ! LIM thermodynamics
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( kideb, kiut )
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 , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
70      !!
71      INTEGER  ::   ji , jk        ! dummy loop indices
72      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
73      INTEGER  ::   iter
74
75      REAL(wp) ::   ztmelts             ! local scalar
76      REAL(wp) ::   zdum       
77      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
78      REAL(wp) ::   zs_snic      ! snow-ice salinity
79      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
80      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
81      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
82      REAL(wp) ::   zgrr         ! bottom growth rate
83      REAL(wp) ::   zt_i_new     ! bottom formation temperature
84
85      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
86      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
87      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
88      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
89      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
90      REAL(wp) ::   zsstK        ! SST in Kelvin
91
92      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3)
93      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2)
94      REAL(wp), POINTER, DIMENSION(:) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
95      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
96      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
97      REAL(wp), POINTER, DIMENSION(:) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
98
99      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt
100      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation
101      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation
102
103      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah
104      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness
105      INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting
106
107      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
108      REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing
109
110      REAL(wp) :: zswitch_sal
111
112      ! Heat conservation
113      INTEGER  ::   num_iter_max
114
115      !!------------------------------------------------------------------
116
117      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
118      SELECT CASE( nn_icesal )                  ! varying salinity or not
119         CASE( 1, 3 ) ;   zswitch_sal = 0       ! prescribed salinity profile
120         CASE( 2 )    ;   zswitch_sal = 1       ! varying salinity profile
121      END SELECT
122
123      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
124      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
125      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
126      CALL wrk_alloc( jpij, nlay_i, icount )
127       
128      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp ; dh_i_sub(:) = 0._wp
129      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
130
131      zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp
132      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp ; zevap_rema(:) = 0._wp ;
133      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp
134
135      zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp       
136      icount (:,:) = 0
137
138
139      ! Initialize enthalpy at nlay_i+1
140      DO ji = kideb, kiut
141         q_i_1d(ji,nlay_i+1) = 0._wp
142      END DO
143
144      ! initialize layer thicknesses and enthalpies
145      h_i_old (:,0:nlay_i+1) = 0._wp
146      qh_i_old(:,0:nlay_i+1) = 0._wp
147      DO jk = 1, nlay_i
148         DO ji = kideb, kiut
149            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
150            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)
151         ENDDO
152      ENDDO
153      !
154      !------------------------------------------------------------------------------!
155      !  1) Calculate available heat for surface and bottom ablation                 !
156      !------------------------------------------------------------------------------!
157      !
158      DO ji = kideb, kiut
159         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
160         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
161
162         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
163         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
164      END DO
165
166      !
167      !------------------------------------------------------------------------------!
168      ! If snow temperature is above freezing point, then snow melts
169      ! (should not happen but sometimes it does)
170      !------------------------------------------------------------------------------!
171      DO ji = kideb, kiut
172         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
173            ! Contribution to heat flux to the ocean [W.m-2], < 0 
174            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
175            ! Contribution to mass flux
176            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
177            ! updates
178            ht_s_1d(ji)   = 0._wp
179            q_s_1d (ji,1) = 0._wp
180            t_s_1d (ji,1) = rt0
181         END IF
182      END DO
183
184      !------------------------------------------------------------!
185      !  2) Computing layer thicknesses and enthalpies.            !
186      !------------------------------------------------------------!
187      !
188      DO jk = 1, nlay_i
189         DO ji = kideb, kiut
190            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
191            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
192         END DO
193      END DO
194      !
195      !------------------------------------------------------------------------------|
196      !  3) Surface ablation and sublimation                                         |
197      !------------------------------------------------------------------------------|
198      !
199      !-------------------------
200      ! 3.1 Snow precips / melt
201      !-------------------------
202      ! Snow accumulation in one thermodynamic time step
203      ! snowfall is partitionned between leads and ice
204      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
205      ! but because of the winds, more snow falls on leads than on sea ice
206      ! and a greater fraction (1-at_i)^beta of the total mass of snow
207      ! (beta < 1) falls in leads.
208      ! In reality, beta depends on wind speed,
209      ! and should decrease with increasing wind speed but here, it is
210      ! considered as a constant. an average value is 0.66
211      ! Martin Vancoppenolle, December 2006
212
213      CALL lim_thd_snwblow( 1. - at_i_1d(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing
214
215      zdeltah(:,:) = 0._wp
216      DO ji = kideb, kiut
217         !-----------
218         ! Snow fall
219         !-----------
220         ! thickness change
221         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
222         ! enthalpy of the precip (>0, J.m-3)
223         zqprec   (ji) = - qprec_ice_1d(ji)   
224         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
225         ! heat flux from snow precip (>0, W.m-2)
226         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
227         ! mass flux, <0
228         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
229
230         !---------------------
231         ! Melt of falling snow
232         !---------------------
233         ! thickness change
234         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
235         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
236         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
237         ! heat used to melt snow (W.m-2, >0)
238         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
239         ! snow melting only = water into the ocean (then without snow precip), >0
240         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   
241         ! updates available heat + precipitations after melting
242         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )     
243         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
244
245         ! update thickness
246         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
247      END DO
248
249      ! If heat still available (zq_su > 0), then melt more snow
250      zdeltah(:,:) = 0._wp
251      DO jk = 1, nlay_s
252         DO ji = kideb, kiut
253            ! thickness change
254            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
255            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 
256            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
257            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting
258            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)   
259            ! heat used to melt snow(W.m-2, >0)
260            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
261            ! snow melting only = water into the ocean (then without snow precip)
262            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
263            ! updates available heat + thickness
264            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
265            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
266         END DO
267      END DO
268
269      !------------------------------
270      ! 3.2 Sublimation (part1: snow)
271      !------------------------------
272      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
273      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean)
274      zdeltah(:,:) = 0._wp
275      DO ji = kideb, kiut
276         zdh_s_sub(ji)  = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
277         ! remaining evap in kg.m-2 (used for ice melting later on)
278         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn
279         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
280         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
281         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  &
282            &                              ) * a_i_1d(ji) * r1_rdtice
283         ! Mass flux by sublimation
284         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
285         ! new snow thickness
286         ht_s_1d(ji)    =  MAX( 0._wp , ht_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 = kideb, kiut
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 = kideb,kiut
303            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
304            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
305              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
306              &              ( ht_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(:,:) = 0._wp ! important
314      DO jk = 1, nlay_i
315         DO ji = kideb, kiut
316            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
317           
318            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
319
320               zEi            = - q_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 sm_i_1d and not s_i_1d(jk) is ok)
335               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_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            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
343               zEw            =    rcp * ( ztmelts - rt0 )            ! 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 sm_i_1d and not s_i_1d(jk) is ok)
361               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_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 * sm_i_1d(ji) * r1_rdtice
382            ! Heat flux [W.m-2], < 0
383            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
384            ! Mass flux > 0
385            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
386            ! update remaining mass flux
387            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
388           
389            ! record which layers have disappeared (for bottom melting)
390            !    => icount=0 : no layer has vanished
391            !    => icount=5 : 5 layers have vanished
392            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
393            icount(ji,jk) = NINT( rswitch )
394            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
395                       
396            ! update heat content (J.m-2) and layer thickness
397            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
398            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
399         END DO
400      END DO
401      ! update ice thickness
402      DO ji = kideb, kiut
403         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
404      END DO
405
406      ! remaining "potential" evap is sent to ocean
407      DO ji = kideb, kiut
408         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
409         wfx_err_sub(ii,ij) = wfx_err_sub(ii,ij) - 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
433
434      ! Iterative procedure
435      DO iter = 1, num_iter_max
436         DO ji = kideb, kiut
437            IF(  zf_tt(ji) < 0._wp  ) THEN
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               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
451
452               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
453                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
454               ! New ice growth
455               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
456
457               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
458               
459               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
460                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
461                  &               + rcp  * ( ztmelts-rt0 )         
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_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
468
469               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
470               
471            ENDIF
472         END DO
473      END DO
474
475      ! Contribution to Energy and Salt Fluxes
476      DO ji = kideb, kiut
477         IF(  zf_tt(ji) < 0._wp  ) THEN
478            ! New ice growth
479                                   
480            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
481
482            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
483           
484            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
485           
486            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
487               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
488               &               + rcp  * ( ztmelts-rt0 )         
489           
490            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
491           
492            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
493           
494            ! Contribution to heat flux to the ocean [W.m-2], >0 
495            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
496
497            ! Total heat flux used in this process [W.m-2], <0 
498            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
499           
500            ! Contribution to salt flux, <0
501            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
502
503            ! Contribution to mass flux, <0
504            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
505
506            ! update heat content (J.m-2) and layer thickness
507            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1)
508            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
509         ENDIF
510      END DO
511
512      !----------------
513      ! 4.2 Basal melt
514      !----------------
515      zdeltah(:,:) = 0._wp ! important
516      DO jk = nlay_i, 1, -1
517         DO ji = kideb, kiut
518            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
519
520               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
521
522               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
523
524                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
525                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
526                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
527                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
528                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
529
530                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
531
532                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
533
534                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
535                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
536
537                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
538                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
539                                   
540                  ! Contribution to mass flux
541                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
542
543                  ! update heat content (J.m-2) and layer thickness
544                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
545                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
546
547               ELSE                               !!! Basal melting
548
549                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
550                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
551                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
552
553                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
554
555                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
556
557                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
558                 
559                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
560
561                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
562
563                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
564
565                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
566
567                  ! Contribution to heat flux to the ocean [W.m-2], <0 
568                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
569
570                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
571                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
572                 
573                  ! Total heat flux used in this process [W.m-2], >0 
574                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
575                 
576                  ! Contribution to mass flux
577                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
578
579                  ! update heat content (J.m-2) and layer thickness
580                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
581                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
582               ENDIF
583           
584            ENDIF
585         END DO
586      END DO
587
588      !-------------------------------------------
589      ! Update temperature, energy
590      !-------------------------------------------
591      DO ji = kideb, kiut
592         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
593      END DO 
594
595      !-------------------------------------------
596      ! 5. What to do with remaining energy
597      !-------------------------------------------
598      ! If heat still available for melting and snow remains, then melt more snow
599      !-------------------------------------------
600      zdeltah(:,:) = 0._wp ! important
601      DO ji = kideb, kiut
602         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
603         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
604         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )
605         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )
606         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
607         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
608         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
609       
610         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2)
611         ! heat used to melt snow
612         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
613         ! Contribution to mass flux
614         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
615         !   
616         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
617         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
618         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
619
620         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
621      END DO
622     
623      !
624      !------------------------------------------------------------------------------|
625      !  6) Snow-Ice formation                                                       |
626      !------------------------------------------------------------------------------|
627      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
628      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
629      DO ji = kideb, kiut
630         !
631         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
632
633         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji)
634         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji)
635
636         ! Salinity of snow ice
637         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
638         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)
639
640         ! entrapment during snow ice formation
641         ! new salinity difference stored (to be used in limthd_sal.F90)
642         IF (  nn_icesal == 2  ) THEN
643            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
644            ! salinity dif due to snow-ice formation
645            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch     
646            ! salinity dif due to bottom growth
647            IF (  zf_tt(ji)  < 0._wp ) THEN
648               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch
649            ENDIF
650         ENDIF
651
652         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
653         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp )    ! <0
654         zsstK          = sst_m(ii,ij) + rt0                               
655         zEw            = rcp * ( zsstK - rt0 )
656         zQm            = zfmdt * zEw 
657         
658         ! Contribution to heat flux
659         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
660
661         ! Contribution to salt flux
662         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
663
664         ! virtual salt flux to keep salinity constant
665         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
666            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_m(ii,ij) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m into the ocean
667               &                            - sm_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  sm_i  from the ocean
668         ENDIF
669         
670         ! Contribution to mass flux
671         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
672         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
673         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
674
675         ! update heat content (J.m-2) and layer thickness
676         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw
677         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
678         
679      END DO
680
681      !
682      !-------------------------------------------
683      ! Update temperature, energy
684      !-------------------------------------------
685      DO ji = kideb, kiut
686         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
687         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
688      END DO
689
690      DO jk = 1, nlay_s
691         DO ji = kideb,kiut
692            ! mask enthalpy
693            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
694            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)
695            ! recalculate t_s_1d from q_s_1d
696            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
697         END DO
698      END DO
699
700      ! --- ensure that a_i = 0 where ht_i = 0 ---
701      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
702     
703      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
704      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
705      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
706      CALL wrk_dealloc( jpij, nlay_i, icount )
707      !
708      !
709   END SUBROUTINE lim_thd_dh
710
711
712   !!--------------------------------------------------------------------------
713   !! INTERFACE lim_thd_snwblow
714   !! ** Purpose :   Compute distribution of precip over the ice
715   !!--------------------------------------------------------------------------
716   SUBROUTINE lim_thd_snwblow_2d( pin, pout )
717      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) )
718      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
719      pout = ( 1._wp - ( pin )**rn_betas )
720   END SUBROUTINE lim_thd_snwblow_2d
721
722   SUBROUTINE lim_thd_snwblow_1d( pin, pout )
723      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
724      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
725      pout = ( 1._wp - ( pin )**rn_betas )
726   END SUBROUTINE lim_thd_snwblow_1d
727
728   
729#else
730   !!----------------------------------------------------------------------
731   !!   Default option                               NO  LIM3 sea-ice model
732   !!----------------------------------------------------------------------
733CONTAINS
734   SUBROUTINE lim_thd_dh          ! Empty routine
735   END SUBROUTINE lim_thd_dh
736#endif
737
738   !!======================================================================
739END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.