source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90 @ 12475

Last change on this file since 12475 was 12475, checked in by dancopsey, 12 months ago
  • Add more print statements.
  • Move away from using snow to ice diagnostics and use a new snow to pond one instead.
File size: 40.6 KB
Line 
1MODULE icethd_dh
2   !!======================================================================
3   !!                       ***  MODULE icethd_dh ***
4   !!   seaice : thermodynamic growth and melt
5   !!======================================================================
6   !! History :       !  2003-05  (M. Vancoppenolle) Original code in 1D
7   !!                 !  2005-06  (M. Vancoppenolle) 3D version
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_thd_dh        : vertical sea-ice growth and melt
15   !!   ice_thd_snwblow   : distribute snow fall between ice and ocean
16  !!----------------------------------------------------------------------
17   USE dom_oce        ! ocean space and time domain
18   USE phycst         ! physical constants
19   USE ice            ! sea-ice: variables
20   USE ice1D          ! sea-ice: thermodynamics variables
21   USE icethd_sal     ! sea-ice: salinity profiles
22   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
26   USE icethd_pnd, ONLY : a_pnd_avail
27   
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   ice_thd_dh        ! called by ice_thd
32   PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbccpl and here
33
34   INTERFACE ice_thd_snwblow
35      MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d
36   END INTERFACE
37
38   !!----------------------------------------------------------------------
39   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
40   !! $Id$
41   !! Software governed by the CeCILL license (see ./LICENSE)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE ice_thd_dh
46      !!------------------------------------------------------------------
47      !!                ***  ROUTINE ice_thd_dh  ***
48      !!
49      !! ** Purpose :   compute ice and snow thickness changes due to growth/melting
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      !!                - Compute available flux of heat for surface ablation
58      !!                - Compute snow and sea ice enthalpies
59      !!                - Surface ablation and sublimation
60      !!                - Bottom accretion/ablation
61      !!                - 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         ! local integer
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      REAL(wp) ::   z1_rho       ! 1/(rhos+rau0-rhoi)
80
81      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
82      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
83      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
84      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
85      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
86
87      REAL(wp), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3)
88      REAL(wp), DIMENSION(jpij) ::   zq_top      ! heat for surface ablation                   (J.m-2)
89      REAL(wp), DIMENSION(jpij) ::   zq_bot      ! heat for bottom ablation                    (J.m-2)
90      REAL(wp), DIMENSION(jpij) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
91      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
92      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
93
94      REAL(wp), DIMENSION(jpij) ::   zdh_s_mel   ! snow melt
95      REAL(wp), DIMENSION(jpij) ::   zdh_s_pre   ! snow precipitation
96      REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation
97
98      REAL(wp), DIMENSION(jpij,nlay_s) ::   zh_s      ! snw layer thickness
99      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
100      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
101      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah_pnd  ! Change in thickness due to melting into ponds
102      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting
103
104      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
105
106      REAL(wp) ::   zswitch_sal
107
108      INTEGER  ::   num_iter_max      ! Heat conservation
109      !!------------------------------------------------------------------
110
111      ! Discriminate between time varying salinity and constant
112      SELECT CASE( nn_icesal )                  ! varying salinity or not
113         CASE( 1 , 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
114         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
115      END SELECT
116
117      ! Initialise fraction of sea ice available for melt ponding
118      DO ji = 1, npti
119         a_pnd_avail_1d(ji) = a_pnd_avail * at_i_1d(ji)
120      END DO
121
122      ! initialize layer thicknesses and enthalpies
123      h_i_old (1:npti,0:nlay_i+1) = 0._wp
124      eh_i_old(1:npti,0:nlay_i+1) = 0._wp
125      DO jk = 1, nlay_i
126         DO ji = 1, npti
127            h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
128            eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
129         END DO
130      END DO
131      !
132      !                       ! ============================================== !
133      !                       ! Available heat for surface and bottom ablation !
134      !                       ! ============================================== !
135      !
136      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
137         !
138         DO ji = 1, npti
139            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
140
141            IF ( to_print(ji) == 10 ) THEN
142               write(numout,*)'loading flux: zq_top(ji), qml_ice_1d(ji) = ',zq_top(ji), ' ', qml_ice_1d(ji)
143            ENDIF
144         END DO
145         !
146      ELSE
147         !
148         DO ji = 1, npti
149            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
150            qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
151            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
152
153            IF ( to_print(ji) == 10 ) THEN
154               write(numout,*)'icethd_dh: qns_ice_1d(ji), qsr_ice_1d(ji), qtr_ice_top_1d(ji), qcn_ice_top_1d(ji) = ',qns_ice_1d(ji), ' ', qsr_ice_1d(ji), ' ', qtr_ice_top_1d(ji), ' ', qcn_ice_top_1d(ji)
155               write(numout,*)'icethd_dh: zdum, t_su_1d(ji), qml_ice_1d(ji), zq_top(ji) = ',zdum, ' ', t_su_1d(ji), qml_ice_1d(ji), zq_top(ji)
156            ENDIF
157
158         END DO
159         !
160      ENDIF
161      !
162      DO ji = 1, npti
163         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 
164         zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rdt_ice )
165         IF ( to_print(ji) == 1 ) THEN
166            write(numout,*)'icethd_dh: qcn_ice_bot_1d(ji), qsb_ice_bot_1d(ji) = ',qcn_ice_bot_1d(ji), qsb_ice_bot_1d(ji)
167            write(numout,*)'icethd_dh: fhld_1d(ji), zq_bot(ji) = ',fhld_1d(ji), ', ', zq_bot(ji)
168         ENDIF
169      END DO
170
171      ! Ice and snow layer thicknesses
172      !-------------------------------
173      DO jk = 1, nlay_i
174         DO ji = 1, npti
175            zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
176         END DO
177      END DO
178
179      DO jk = 1, nlay_s
180         DO ji = 1, npti
181            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
182         END DO
183      END DO
184
185      !                       ! ============ !
186      !                       !     Snow     !
187      !                       ! ============ !
188      !
189      ! Internal melting
190      ! ----------------
191      ! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
192      DO jk = 1, nlay_s
193         DO ji = 1, npti
194            IF( t_s_1d(ji,jk) > rt0 ) THEN
195               hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! heat flux to the ocean [W.m-2], < 0
196               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux
197               ! updates
198               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk)
199               h_s_1d  (ji)    = h_s_1d(ji) - zh_s(ji,jk)
200               zh_s    (ji,jk) = 0._wp
201               e_s_1d  (ji,jk) = 0._wp
202               t_s_1d  (ji,jk) = rt0
203            END IF
204
205           IF (to_print(ji) == 10) THEN
206                write(numout,*)'icethd_dh 8: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
207            END IF
208         END DO
209      END DO
210
211      DO ji = 1, npti
212         IF (to_print(ji) == 10) THEN
213            write(numout,*)'icethd_dh internal_melting: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
214         END IF
215      END DO       
216
217      ! Snow precipitation
218      !-------------------
219      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
220
221      zdeltah(1:npti,:) = 0._wp
222      DO ji = 1, npti
223         IF( sprecip_1d(ji) > 0._wp ) THEN
224            !
225            ! --- precipitation ---
226            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)   ! thickness change
227            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3)
228            !
229            hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat flux from snow precip (>0, W.m-2)
230            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0
231           
232            ! --- melt of falling snow ---
233            rswitch              = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
234            zdeltah       (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change
235            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                 ! bound melting
236            hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat used to melt snow (W.m-2, >0)
237            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip), >0
238           
239            ! updates available heat + precipitations after melting
240            dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)
241            zq_top   (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )     
242            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
243
244            IF ( to_print(ji) == 10 ) THEN
245               write(numout,*)'after precip: zq_top(ji), zdeltah(ji,1), zqprec(ji) = ',zq_top(ji), zdeltah(ji,1), zqprec(ji)
246            ENDIF
247           
248            ! update thickness
249            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
250            !
251         ELSE
252            !
253            zdh_s_pre(ji) = 0._wp
254            zqprec   (ji) = 0._wp
255            !
256         ENDIF
257
258         IF (to_print(ji) == 10) THEN
259            write(numout,*)'icethd_dh snow_precip: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
260         END IF
261
262      END DO
263
264      ! recalculate snow layers
265      DO jk = 1, nlay_s
266         DO ji = 1, npti
267            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
268         END DO
269      END DO
270
271      ! Snow melting
272      ! ------------
273      ! If heat still available (zq_top > 0), then melt more snow
274      zdeltah(1:npti,:) = 0._wp
275      zdh_s_mel(1:npti) = 0._wp
276      DO jk = 1, nlay_s
277         DO ji = 1, npti
278            IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
279               !
280               rswitch          = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) )
281               zdeltah  (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change
282               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting
283               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)
284               zdeltah_pnd(ji,jk) = zdeltah  (ji,jk) * a_pnd_avail_1d(ji)
285               
286               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0)
287               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * (zdeltah(ji,jk) - zdeltah_pnd(ji,jk)) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip)
288               
289               ! updates available heat + thickness
290               dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk)
291               dh_s_pnd(ji)    = dh_s_pnd(ji) + zdeltah_pnd(ji,jk)                ! Cumulate surface melt on areas contributing to melt ponds
292               zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
293               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
294               zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
295
296               ! Snow melt trapped in ponds contrbutes to ice thickness and snow to ice conversion
297               h_i_1d  (ji)    = MAX( 0._wp , h_i_1d(ji) - zdeltah_pnd(ji,jk) * rhos * r1_rhoi ) 
298               wfx_snp_1d(ji)     = wfx_snp_1d(ji)     + a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice  ! As zdeltah_pnd is negative this comes out as negative
299               !
300            ENDIF
301
302            IF ( to_print(ji) == 10 ) THEN
303               write(numout,*)'after snow melt: zq_top(ji), zdeltah(ji,jk), e_s_1d(ji,jk), test = ',zq_top(ji), zdeltah(ji,jk), e_s_1d(ji,jk), ( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp )
304            ENDIF
305
306            IF (to_print(ji) == 10) THEN
307               write(numout,*)'icethd_dh loop: wfx_snw_sum_1d(ji), jk = ',wfx_snw_sum_1d(ji), ' ', jk
308            END IF
309         END DO
310      END DO
311
312      DO ji = 1, npti
313         IF (to_print(ji) == 10) THEN
314            write(numout,*)'icethd_dh snow_ground: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
315         END IF
316      END DO
317
318      ! Snow sublimation
319      !-----------------
320      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
321      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
322      zdeltah(1:npti,:) = 0._wp
323      DO ji = 1, npti
324         IF( evap_ice_1d(ji) > 0._wp ) THEN
325            !
326            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice )
327            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos   ! remaining evap in kg.m-2 (used for ice melting later on)
328            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
329           
330            hfx_sub_1d    (ji) = hfx_sub_1d(ji) + &   ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow)
331               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  &
332               &                 * a_i_1d(ji) * r1_rdtice
333            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice   ! Mass flux by sublimation
334           
335            ! new snow thickness
336            h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
337            ! update precipitations after sublimation and correct sublimation
338            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
339            zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
340            !
341         ELSE
342            !
343            zdh_s_sub (ji) = 0._wp
344            zevap_rema(ji) = 0._wp
345            !
346         ENDIF
347      END DO
348     
349      ! --- Update snow diags --- !
350      DO ji = 1, npti
351         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
352      END DO
353
354      ! Update temperature, energy
355      !---------------------------
356      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
357      DO jk = 1, nlay_s
358         DO ji = 1,npti
359            rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
360            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            &
361              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  &
362              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) )
363
364            IF (to_print(ji) == 10) THEN
365                write(numout,*)'icethd_dh 9: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
366                write(numout,*)'icethd_dh 9: zdh_s_pre(ji), zqprec(ji), h_s_1d(ji) = ',zdh_s_pre(ji), ' ', zqprec(ji), ' ', h_s_1d(ji)
367            END IF
368         END DO
369      END DO
370     
371      !                       ! ============ !
372      !                       !     Ice      !
373      !                       ! ============ !
374
375      ! Surface ice melting
376      !--------------------
377      zdeltah(1:npti,:) = 0._wp ! important
378      DO jk = 1, nlay_i
379         DO ji = 1, npti
380            ztmelts = - rTmlt * sz_i_1d(ji,jk)   ! Melting point of layer k [C]
381           
382            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
383
384               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]       
385               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0)
386                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
387               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
388                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
389               zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0
390                         
391               dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting
392               
393               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
394
395               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0
396               !                                                                                                  ice enthalpy zEi is "sent" to the ocean
397               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux
398               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
399               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
400
401            ELSE                                        !-- Surface melting
402               
403               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]
404               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
405               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
406               
407               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean or meltponds [kg/m2, >0]
408               
409               zdeltah(ji,jk) = - zfmdt * r1_rhoi * (1._wp - a_pnd_avail_1d(ji))                    ! Melt of layer jk [m, <0] into the ocean
410               zdeltah_pnd(ji,jk) = - zfmdt * r1_rhoi * a_pnd_avail_1d(ji)                          ! Melt of layer jk [m, <0] into meltponds
411               
412               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]
413               
414               zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat
415               
416               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)          ! Cumulate surface melt on areas not contributing to melt ponds
417               dh_i_pnd(ji)   = dh_i_pnd(ji) + zdeltah_pnd(ji,jk)      ! Cumulate surface melt on areas contributing to melt ponds
418               
419               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
420               
421               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
422               
423               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux >0
424               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok)
425               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0
426               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], > 0 
427               !
428               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
429               
430            END IF
431
432            IF ( to_print(ji) == 10 ) THEN
433               write(numout,*)'after top melt: zq_top(ji), test = ',zq_top(ji), ( t_i_1d(ji,jk) >= (ztmelts+rt0) ) 
434               write(numout,*)'after top melt: wfx_sum_1d(ji), zdeltah(ji,jk), zfmdt, zdE = ', wfx_sum_1d(ji), zdeltah(ji,jk), zfmdt, zdE
435            ENDIF
436           
437            ! Ice sublimation
438            ! ---------------
439            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi )
440            zdeltah (ji,jk) = zdeltah (ji,jk) + zdum
441            dh_i_sub(ji)    = dh_i_sub(ji)    + zdum
442           
443            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice  ! Salt flux >0
444            !                                                                                          clem: flux is sent to the ocean for simplicity
445            !                                                                                                but salt should remain in the ice except
446            !                                                                                                if all ice is melted. => must be corrected
447            hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice      ! Heat flux [W.m-2], < 0
448
449            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_rdtice           ! Mass flux > 0
450
451            ! update remaining mass flux
452            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi
453           
454            ! record which layers have disappeared (for bottom melting)
455            !    => icount=0 : no layer has vanished
456            !    => icount=5 : 5 layers have vanished
457            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
458            icount(ji,jk) = NINT( rswitch )
459            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
460                       
461            ! update heat content (J.m-2) and layer thickness
462            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
463            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
464         END DO
465      END DO
466     
467      ! update ice thickness
468      DO ji = 1, npti
469         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) )
470      END DO
471
472      ! remaining "potential" evap is sent to ocean
473      DO ji = 1, npti
474         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)
475      END DO
476
477
478      ! Ice Basal growth
479      !------------------
480      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
481      ! between the inner conductive flux  (qcn_ice_bot), from the open water heat flux
482      ! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
483      ! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice
484
485      ! If salinity varies in time, an iterative procedure is required, because
486      ! the involved quantities are inter-dependent.
487      ! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
488      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
489      ! -> need for an iterative procedure, which converges quickly
490
491      num_iter_max = 1
492      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
493     
494      DO ji = 1, npti
495         IF(  zf_tt(ji) < 0._wp  ) THEN
496            DO iter = 1, num_iter_max   ! iterations
497
498               ! New bottom ice salinity (Cox & Weeks, JGR88 )
499               !--- zswi1  if dh/dt < 2.0e-8
500               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
501               !--- zswi2  if dh/dt > 3.6e-7
502               zgrr     = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_rdtice , epsi10 ) )
503               zswi2    = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
504               zswi12   = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
505               zswi1    = 1. - zswi2 * zswi12
506               zfracs   = MIN( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
507                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
508
509               s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity
510
511               ztmelts       = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C)
512
513               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
514               
515               zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
516                  &            - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts
517
518               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
519
520               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0)
521
522               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
523               
524            END DO
525            ! Contribution to Energy and Salt Fluxes                                   
526            zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0)
527           
528            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], >0
529            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], <0
530           
531            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice     ! Salt flux, <0
532
533            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice                   ! Mass flux, <0
534
535            ! update heat content (J.m-2) and layer thickness
536            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
537            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
538
539         ENDIF
540
541      END DO
542
543      ! Ice Basal melt
544      !---------------
545      zdeltah(1:npti,:) = 0._wp ! important
546      DO jk = nlay_i, 1, -1
547         DO ji = 1, npti
548            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
549
550               ztmelts = - rTmlt * sz_i_1d(ji,jk)  ! Melting point of layer jk (C)
551
552               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
553
554                  zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0)
555                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
556                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted)
557                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
558                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
559
560                  dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk)
561
562                  zfmdt             = - zdeltah(ji,jk) * rhoi      ! Mass flux x time step > 0
563
564                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0
565                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean
566                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux
567                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
568                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
569
570                  ! update heat content (J.m-2) and layer thickness
571                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
572                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
573
574               ELSE                                        !-- Basal melting
575
576                  zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0)
577                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0)
578                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0)
579
580                  zfmdt           = - zq_bot(ji) / zdE                                        ! Mass flux x time step (kg/m2, >0)
581
582                  zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change
583
584                  IF ( to_print(ji) == 1 ) THEN
585                    write(numout,*)'icethd_dh: zfmdt, zq_bot(ji), zdE, zdeltah(ji,jk) = ',zfmdt,', ', zq_bot(ji), zdE, zdeltah(ji,jk)
586                  ENDIF
587
588                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change
589                 
590                  zq_bot(ji)      = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors
591
592                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt
593
594                  zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0
595
596                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean
597
598                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
599                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat used in this process [W.m-2], >0 
600
601                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux
602                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok
603                 
604                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
605
606                  IF (to_print(ji) == 1) THEN
607                     write(numout,*)'icethd_dh: wfx_bom_1d(ji), in meters, at_i_1d(ji) = ',wfx_bom_1d(ji),wfx_bom_1d(ji)*rdt_ice/(rhoi*at_i_1d(ji)), at_i_1d(ji)
608                  ENDIF
609
610                  ! update heat content (J.m-2) and layer thickness
611                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
612                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
613               ENDIF
614           
615            ENDIF
616         END DO
617      END DO
618
619     
620
621      ! Update temperature, energy
622      ! --------------------------
623      DO ji = 1, npti
624         h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) )
625      END DO 
626
627      ! If heat still available then melt more snow
628      !-------------------------------------------
629      zdeltah(1:npti,:) = 0._wp ! important
630      DO ji = 1, npti
631         zq_rema (ji)   = zq_top(ji) + zq_bot(ji) 
632         rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
633         rswitch        = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
634         zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
635         zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
636         dh_s_tot(ji)   = dh_s_tot(ji) + zdeltah(ji,1)
637         h_s_1d  (ji)   = h_s_1d  (ji) + zdeltah(ji,1)
638       
639         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2)
640         hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice   ! Heat used to melt snow, W.m-2 (>0)
641         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice       ! Mass flux
642         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
643         !   
644         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
645         qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
646
647         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
648
649         IF (to_print(ji) == 10) THEN
650            write(numout,*)'icethd_dh snow_extra: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
651         END IF
652
653      END DO
654
655      !
656      ! Snow-Ice formation
657      ! ------------------
658      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
659      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
660      z1_rho = 1._wp / ( rhos+rau0-rhoi )
661      DO ji = 1, npti
662         !
663         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho )
664
665         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
666         h_s_1d(ji)    = h_s_1d(ji) - dh_snowice(ji)
667
668         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
669         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0
670         zEw            = rcp * sst_1d(ji)
671         zQm            = zfmdt * zEw 
672         
673         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux
674
675         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux
676
677         ! Case constant salinity in time: virtual salt flux to keep salinity constant
678         IF( nn_icesal /= 2 )  THEN
679            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
680               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice     ! and get  rn_icesal from the ocean
681         ENDIF
682
683         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
684         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice
685         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_rdtice
686
687         ! update heat content (J.m-2) and layer thickness
688         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
689         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
690         
691      END DO
692
693      !
694      ! Update temperature, energy
695      ! --------------------------
696      DO ji = 1, npti
697         rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
698         t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
699
700         IF (to_print(ji) == 10) THEN
701            write(numout,*)'icethd_dh: t_su_1d(ji), t_i_1d(ji,1) = ',t_su_1d(ji), ' ', t_i_1d(ji,1)
702         END IF
703      END DO
704
705      DO jk = 1, nlay_s
706         DO ji = 1,npti
707            ! where there is no ice or no snow
708            rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) )
709            ! mass & energy loss to the ocean
710            hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * &
711               &                              ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! heat flux to the ocean [W.m-2], < 0
712            wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * &
713               &                              ( rhos          * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! mass flux
714            ! update energy (mass is updated in the next loop)
715            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
716            ! recalculate t_s_1d from e_s_1d
717            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
718
719            IF (to_print(ji) == 10) THEN
720                write(numout,*)'icethd_dh 10: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
721            END IF
722         END DO
723      END DO
724
725      ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
726      WHERE( h_i_1d(1:npti) == 0._wp )   
727         a_i_1d(1:npti) = 0._wp
728         h_s_1d(1:npti) = 0._wp
729      END WHERE
730      !
731   END SUBROUTINE ice_thd_dh
732
733
734   !!--------------------------------------------------------------------------
735   !! INTERFACE ice_thd_snwblow
736   !!
737   !! ** Purpose :   Compute distribution of precip over the ice
738   !!
739   !!                Snow accumulation in one thermodynamic time step
740   !!                snowfall is partitionned between leads and ice.
741   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
742   !!                but because of the winds, more snow falls on leads than on sea ice
743   !!                and a greater fraction (1-at_i)^beta of the total mass of snow
744   !!                (beta < 1) falls in leads.
745   !!                In reality, beta depends on wind speed,
746   !!                and should decrease with increasing wind speed but here, it is
747   !!                considered as a constant. an average value is 0.66
748   !!--------------------------------------------------------------------------
749!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
750   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
751      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
752      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
753      pout = ( 1._wp - ( pin )**rn_blow_s )
754   END SUBROUTINE ice_thd_snwblow_2d
755
756   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
757      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
758      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
759      pout = ( 1._wp - ( pin )**rn_blow_s )
760   END SUBROUTINE ice_thd_snwblow_1d
761
762#else
763   !!----------------------------------------------------------------------
764   !!   Default option                                NO SI3 sea-ice model
765   !!----------------------------------------------------------------------
766#endif
767
768   !!======================================================================
769END MODULE icethd_dh
Note: See TracBrowser for help on using the repository browser.