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.
limthd_dh.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 8369 was 8369, checked in by clem, 5 years ago

STEP4 (4): put all thermodynamics in 1D (limitd_th OK)

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