source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 4 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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