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.
icethd_dh.F90 in NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE – NEMO

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

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