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_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 36.3 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 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( kideb, kiut )
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 , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
69      !!
70      INTEGER  ::   ji , jk        ! dummy loop indices
71      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
72      INTEGER  ::   iter
73
74      REAL(wp) ::   ztmelts             ! local scalar
75      REAL(wp) ::   zdum       
76      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
77      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
78      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
79      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
80      REAL(wp) ::   zgrr         ! bottom growth rate
81      REAL(wp) ::   zt_i_new     ! bottom formation temperature
82
83      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
84      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
85      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
86      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
87      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
88      REAL(wp) ::   zsstK        ! SST in Kelvin
89
90      REAL(wp), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3)
91      REAL(wp), DIMENSION(jpij) ::   zq_su       ! heat for surface ablation                   (J.m-2)
92      REAL(wp), DIMENSION(jpij) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
93      REAL(wp), DIMENSION(jpij) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
94      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
95      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
96
97      REAL(wp), DIMENSION(jpij) ::   zdh_s_mel   ! snow melt
98      REAL(wp), DIMENSION(jpij) ::   zdh_s_pre   ! snow precipitation
99      REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation
100
101      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
102      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
103      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting
104
105      REAL(wp), DIMENSION(jpij) ::   zqh_i       ! total ice heat content  (J.m-2)
106      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
107
108      REAL(wp) :: zswitch_sal
109
110      ! Heat conservation
111      INTEGER  ::   num_iter_max
112
113      !!------------------------------------------------------------------
114
115      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
116      SELECT CASE( nn_icesal )                  ! varying salinity or not
117         CASE( 1, 3 ) ;   zswitch_sal = 0       ! prescribed salinity profile
118         CASE( 2 )    ;   zswitch_sal = 1       ! varying salinity profile
119      END SELECT
120
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 ; zqh_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 = kideb, kiut
131         q_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      qh_i_old(:,0:nlay_i+1) = 0._wp
137      DO jk = 1, nlay_i
138         DO ji = kideb, kiut
139            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
140            qh_i_old(ji,jk) = q_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 = kideb, kiut
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 = kideb, kiut
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) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
165            ! Contribution to mass flux
166            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
167            ! updates
168            ht_s_1d(ji)   = 0._wp
169            q_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 = kideb, kiut
180            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
181            zqh_i(ji)   = zqh_i(ji) + q_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(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing
204
205      zdeltah(:,:) = 0._wp
206      DO ji = kideb, kiut
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_1d(ji) = wfx_snw_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 = kideb, kiut
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, q_s_1d(ji,jk) - epsi20 ) ) ) 
246            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_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) * q_s_1d(ji,jk) * r1_rdtice 
251            ! snow melting only = water into the ocean (then without snow precip)
252            wfx_snw_1d(ji)   = wfx_snw_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) * q_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 = kideb, kiut
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) ) * q_s_1d(ji,1)  &
272            &                              ) * a_i_1d(ji) * r1_rdtice
273         ! Mass flux by sublimation
274         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
275         ! new snow thickness
276         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
277         ! update precipitations after sublimation and correct sublimation
278         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
279         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
280      END DO
281     
282      ! --- Update snow diags --- !
283      DO ji = kideb, kiut
284         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
285      END DO
286
287      !-------------------------------------------
288      ! 3.3 Update temperature, energy
289      !-------------------------------------------
290      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
291      DO jk = 1, nlay_s
292         DO ji = kideb,kiut
293            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
294            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
295              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
296              &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
297         END DO
298      END DO
299
300      !--------------------------
301      ! 3.4 Surface ice ablation
302      !--------------------------
303      zdeltah(:,:) = 0._wp ! important
304      DO jk = 1, nlay_i
305         DO ji = kideb, kiut
306            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
307           
308            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
309
310               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
311               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
312                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
313               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
314                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
315               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0
316                         
317               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
318               
319               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
320
321               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
322               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
323               
324               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
325               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
326               
327               ! Contribution to mass flux
328               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
329
330            ELSE                                !!! Surface melting
331               
332               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
333               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
334               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
335               
336               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
337               
338               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
339               
340               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]
341               
342               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
343               
344               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
345               
346               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
347               
348               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
349               
350               ! Contribution to salt flux >0 (clem: using sm_i_1d and not s_i_1d(jk) is ok)
351               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
352               
353               ! Contribution to heat flux [W.m-2], < 0
354               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
355               
356               ! Total heat flux used in this process [W.m-2], > 0 
357               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
358               
359               ! Contribution to mass flux
360               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
361               
362            END IF
363            ! ----------------------
364            ! Sublimation part2: ice
365            ! ----------------------
366            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
367            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum
368            dh_i_sub(ji)  = dh_i_sub(ji) + zdum
369            ! 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.
370            !                          It must be corrected at some point)
371            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice
372            ! Heat flux [W.m-2], < 0
373            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
374            ! Mass flux > 0
375            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
376            ! update remaining mass flux
377            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
378           
379            ! record which layers have disappeared (for bottom melting)
380            !    => icount=0 : no layer has vanished
381            !    => icount=5 : 5 layers have vanished
382            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
383            icount(ji,jk) = NINT( rswitch )
384            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
385                       
386            ! update heat content (J.m-2) and layer thickness
387            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
388            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
389         END DO
390      END DO
391      ! update ice thickness
392      DO ji = kideb, kiut
393         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
394      END DO
395
396      ! remaining "potential" evap is sent to ocean
397      DO ji = kideb, kiut
398         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
399         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)
400      END DO
401
402      !
403      !------------------------------------------------------------------------------!
404      ! 4) Basal growth / melt                                                       !
405      !------------------------------------------------------------------------------!
406      !
407      !------------------
408      ! 4.1 Basal growth
409      !------------------
410      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
411      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux
412      ! (fhld) and the turbulent ocean flux (fhtur).
413      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
414
415      ! If salinity varies in time, an iterative procedure is required, because
416      ! the involved quantities are inter-dependent.
417      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
418      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
419      ! -> need for an iterative procedure, which converges quickly
420
421      num_iter_max = 1
422      IF( nn_icesal == 2 ) num_iter_max = 5
423
424      ! Iterative procedure
425      DO iter = 1, num_iter_max
426         DO ji = kideb, kiut
427            IF(  zf_tt(ji) < 0._wp  ) THEN
428
429               ! New bottom ice salinity (Cox & Weeks, JGR88 )
430               !--- zswi1  if dh/dt < 2.0e-8
431               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
432               !--- zswi2  if dh/dt > 3.6e-7
433               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
434               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
435               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
436               zswi1              = 1. - zswi2 * zswi12
437               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
438                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
439
440               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
441
442               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! 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               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
460               
461            ENDIF
462         END DO
463      END DO
464
465      ! Contribution to Energy and Salt Fluxes
466      DO ji = kideb, kiut
467         IF(  zf_tt(ji) < 0._wp  ) THEN
468            ! New ice growth
469                                   
470            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
471
472            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
473           
474            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
475           
476            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
477               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
478               &               + rcp  * ( ztmelts-rt0 )         
479           
480            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
481           
482            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
483           
484            ! Contribution to heat flux to the ocean [W.m-2], >0 
485            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
486
487            ! Total heat flux used in this process [W.m-2], <0 
488            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
489           
490            ! Contribution to salt flux, <0
491            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
492
493            ! Contribution to mass flux, <0
494            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
495
496            ! update heat content (J.m-2) and layer thickness
497            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)
498            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
499         ENDIF
500      END DO
501
502      !----------------
503      ! 4.2 Basal melt
504      !----------------
505      zdeltah(:,:) = 0._wp ! important
506      DO jk = nlay_i, 1, -1
507         DO ji = kideb, kiut
508            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
509
510               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
511
512               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
513
514                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
515                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
516                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
517                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
518                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
519
520                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
521
522                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
523
524                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
525                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
526
527                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
528                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
529                                   
530                  ! Contribution to mass flux
531                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
532
533                  ! update heat content (J.m-2) and layer thickness
534                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
535                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
536
537               ELSE                               !!! Basal melting
538
539                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
540                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
541                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
542
543                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
544
545                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
546
547                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
548                 
549                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
550
551                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
552
553                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
554
555                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
556
557                  ! Contribution to heat flux to the ocean [W.m-2], <0 
558                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
559
560                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
561                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
562                 
563                  ! Total heat flux used in this process [W.m-2], >0 
564                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
565                 
566                  ! Contribution to mass flux
567                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
568
569                  ! update heat content (J.m-2) and layer thickness
570                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
571                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
572               ENDIF
573           
574            ENDIF
575         END DO
576      END DO
577
578      !-------------------------------------------
579      ! Update temperature, energy
580      !-------------------------------------------
581      DO ji = kideb, kiut
582         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
583      END DO 
584
585      !-------------------------------------------
586      ! 5. What to do with remaining energy
587      !-------------------------------------------
588      ! If heat still available for melting and snow remains, then melt more snow
589      !-------------------------------------------
590      zdeltah(:,:) = 0._wp ! important
591      DO ji = kideb, kiut
592         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
593         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
594         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )
595         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )
596         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
597         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
598         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
599       
600         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2)
601         ! heat used to melt snow
602         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)
603         ! Contribution to mass flux
604         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
605         !   
606         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
607         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
608         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
609
610         IF( ln_limctl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
611      END DO
612     
613      !
614      !------------------------------------------------------------------------------|
615      !  6) Snow-Ice formation                                                       |
616      !------------------------------------------------------------------------------|
617      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
618      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
619      DO ji = kideb, kiut
620         !
621         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
622
623         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji)
624         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji)
625
626         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
627         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
628         zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0
629         zsstK          = sst_m(ii,ij) + rt0                               
630         zEw            = rcp * ( zsstK - rt0 )
631         zQm            = zfmdt * zEw 
632         
633         ! Contribution to heat flux
634         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
635
636         ! Contribution to salt flux
637         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
638
639         ! virtual salt flux to keep salinity constant
640         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
641            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
642               &                            - sm_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean
643         ENDIF
644
645         ! Contribution to mass flux
646         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
647         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
648         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
649
650         ! update heat content (J.m-2) and layer thickness
651         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw
652         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
653         
654      END DO
655
656      !
657      !-------------------------------------------
658      ! Update temperature, energy
659      !-------------------------------------------
660      DO ji = kideb, kiut
661         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
662         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
663      END DO
664
665      DO jk = 1, nlay_s
666         DO ji = kideb,kiut
667            ! mask enthalpy
668            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
669            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)
670            ! recalculate t_s_1d from q_s_1d
671            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
672         END DO
673      END DO
674
675      ! --- ensure that a_i = 0 where ht_i = 0 ---
676      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
677     
678      !
679      !
680   END SUBROUTINE lim_thd_dh
681
682
683   !!--------------------------------------------------------------------------
684   !! INTERFACE lim_thd_snwblow
685   !! ** Purpose :   Compute distribution of precip over the ice
686   !!--------------------------------------------------------------------------
687   SUBROUTINE lim_thd_snwblow_2d( pin, pout )
688      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) )
689      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
690      pout = ( 1._wp - ( pin )**rn_betas )
691   END SUBROUTINE lim_thd_snwblow_2d
692
693   SUBROUTINE lim_thd_snwblow_1d( pin, pout )
694      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
695      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
696      pout = ( 1._wp - ( pin )**rn_betas )
697   END SUBROUTINE lim_thd_snwblow_1d
698
699   
700#else
701   !!----------------------------------------------------------------------
702   !!   Default option                               NO  LIM3 sea-ice model
703   !!----------------------------------------------------------------------
704CONTAINS
705   SUBROUTINE lim_thd_dh          ! Empty routine
706   END SUBROUTINE lim_thd_dh
707#endif
708
709   !!======================================================================
710END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.