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

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

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

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

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

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

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