New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limthd_dh.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 9270

Last change on this file since 9270 was 7993, checked in by frrh, 7 years ago

Merge in missing revisions 6428:2477 inclusive and 6482 from nemo_v3_6_STABLE
branch. In ptic, this includes the fix for restartability of runoff fields in coupled
models. Evolution of coupled models will therefor be affected.

These changes donot affect evolution of the current stand-alone NEMO-CICE GO6
standard configuration.

Work and testing documented in Met Office GMED ticket 320.

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