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.F90 in NEMO/branches/UKMO/NEMO_4.0.4_remade_heat_balance/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_remade_heat_balance/src/ICE/icethd.F90 @ 15814

Last change on this file since 15814 was 14585, checked in by dancopsey, 3 years ago
  • Calculate the heat flux used up in lateral melting (hfx_lam). This initial incarnation is not quite right as it released all the enthalpy to the ocean and not just the enthalpy required to heat the ice to the SST.
  • Rearrange the calculation of qt_oce_ai in iceupdate.F90. This removes some of the heat fluxes caused by changes that are not accounted for elswhere in the energy balance such as snow precip and sublimation. It also replaces hfx_thd with hfx_lam. Most of the contributions to hfx_thd were for energy adjustments for all the heat flux terms to represent the temperature difference between the SST and 0oC. These are not needed. However hfx_thd also included a contribution for lateral melt which has now been calculated separately and will be included with all the other ocean to sea ice fluxes.
File size: 36.0 KB
Line 
1MODULE icethd
2   !!======================================================================
3   !!                  ***  MODULE icethd   ***
4   !!   sea-ice : master routine for thermodynamics
5   !!======================================================================
6   !! History :  1.0  !  2000-01  (M.A. Morales Maqueda, H. Goosse, T. Fichefet) original code 1D
7   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_thd       : thermodynamics of sea ice
14   !!   ice_thd_init  : initialisation of sea-ice thermodynamics
15   !!----------------------------------------------------------------------
16   USE phycst         ! physical constants
17   USE dom_oce        ! ocean space and time domain variables
18   USE ice            ! sea-ice: variables
19!!gm list trop longue ==>>> why not passage en argument d'appel ?
20   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl
21   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
22      &                 qml_ice, qcn_ice, qtr_ice_top
23   USE ice1D          ! sea-ice: thermodynamics variables
24   USE icethd_zdf     ! sea-ice: vertical heat diffusion
25   USE icethd_dh      ! sea-ice: ice-snow growth and melt
26   USE icethd_da      ! sea-ice: lateral melting
27   USE icethd_sal     ! sea-ice: salinity
28   USE icethd_ent     ! sea-ice: enthalpy redistribution
29   USE icethd_do      ! sea-ice: growth in open water
30   USE icethd_pnd     ! sea-ice: melt ponds
31   USE iceitd         ! sea-ice: remapping thickness distribution
32   USE icecor         ! sea-ice: corrections
33   USE icetab         ! sea-ice: 1D <==> 2D transformation
34   USE icevar         ! sea-ice: operations
35   USE icectl         ! sea-ice: control print
36   !
37   USE in_out_manager ! I/O manager
38   USE iom            ! I/O manager library
39   USE lib_mpp        ! MPP library
40   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
41   USE lbclnk         ! lateral boundary conditions (or mpp links)
42   USE timing         ! Timing
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   ice_thd         ! called by limstp module
48   PUBLIC   ice_thd_init    ! called by ice_init
49
50   !!** namelist (namthd) **
51   LOGICAL ::   ln_icedH         ! activate ice thickness change from growing/melting (T) or not (F)
52   LOGICAL ::   ln_icedA         ! activate lateral melting param. (T) or not (F)
53   LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F)
54   LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
55   LOGICAL ::   ln_leadhfx       ! heat in the leads is used to melt sea-ice before warming the ocean
56
57   !! for convergence tests
58   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztice_cvgerr, ztice_cvgstp
59
60   !! * Substitutions
61#  include "vectopt_loop_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
64   !! $Id$
65   !! Software governed by the CeCILL license (see ./LICENSE)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE ice_thd( kt )
70      !!-------------------------------------------------------------------
71      !!                ***  ROUTINE ice_thd  ***       
72      !! 
73      !! ** Purpose : This routine manages ice thermodynamics
74      !!         
75      !! ** Action : - computation of oceanic sensible heat flux at the ice base
76      !!                              energy budget in the leads
77      !!                              net fluxes on top of ice and of ocean
78      !!             - selection of grid cells with ice
79      !!                - call ice_thd_zdf  for vertical heat diffusion
80      !!                - call ice_thd_dh   for vertical ice growth and melt
81      !!                - call ice_thd_pnd  for melt ponds
82      !!                - call ice_thd_ent  for enthalpy remapping
83      !!                - call ice_thd_sal  for ice desalination
84      !!                - call ice_thd_temp to  retrieve temperature from ice enthalpy
85      !!                - call ice_thd_mono for extra lateral ice melt if active virtual thickness distribution
86      !!                - call ice_thd_da   for lateral ice melt
87      !!             - back to the geographic grid
88      !!                - call ice_thd_rem  for remapping thickness distribution
89      !!                - call ice_thd_do   for ice growth in leads
90      !!-------------------------------------------------------------------
91      INTEGER, INTENT(in) :: kt    ! number of iteration
92      !
93      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
94      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos
95      REAL(wp), PARAMETER :: zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04)
96      REAL(wp), PARAMETER :: zch        = 0.0057_wp   ! heat transfer coefficient
97      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
98      !
99      !!-------------------------------------------------------------------
100      ! controls
101      IF( ln_timing    )   CALL timing_start('icethd')                                                             ! timing
102      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
103      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
104
105      IF( kt == nit000 .AND. lwp ) THEN
106         WRITE(numout,*)
107         WRITE(numout,*) 'ice_thd: sea-ice thermodynamics'
108         WRITE(numout,*) '~~~~~~~'
109      ENDIF
110
111      ! convergence tests
112      IF( ln_zdf_chkcvg ) THEN
113         ALLOCATE( ztice_cvgerr(jpi,jpj,jpl) , ztice_cvgstp(jpi,jpj,jpl) )
114         ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp
115      ENDIF
116     
117      !---------------------------------------------!
118      ! computation of friction velocity at T points
119      !---------------------------------------------!
120      IF( ln_icedyn ) THEN
121         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
122         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
123         DO jj = 2, jpjm1 
124            DO ji = fs_2, fs_jpim1
125               zfric(ji,jj) = rn_cio * ( 0.5_wp *  &
126                  &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
127                  &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
128               zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + &
129                  &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) )
130            END DO
131         END DO
132      ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean
133         DO jj = 2, jpjm1
134            DO ji = fs_2, fs_jpim1
135               zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp *  &
136                  &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
137                  &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
138               zvel(ji,jj) = 0._wp
139            END DO
140         END DO
141      ENDIF
142      CALL lbc_lnk_multi( 'icethd', zfric, 'T',  1._wp, zvel, 'T', 1._wp )
143      !
144      !--------------------------------------------------------------------!
145      ! Partial computation of forcing for the thermodynamic sea ice model
146      !--------------------------------------------------------------------!
147      DO jj = 1, jpj
148         DO ji = 1, jpi
149            rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
150            !
151            ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- !
152            zqld =  tmask(ji,jj,1) * rdt_ice *  &
153               &    ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) +  &
154               &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
155
156            ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- !
157            !     (mostly<0 but >0 if supercooling)
158            zqfr     = rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)
159            zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0
160            zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0
161
162            ! --- Sensible ocean-to-ice heat flux (W/m2) --- !
163            !     (mostly>0 but <0 if supercooling)
164            zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
165            qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )
166
167            ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach
168            !                              the freezing point, so that we do not have SST < T_freeze
169            !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg
170            !                              The following formulation is ok for both normal conditions and supercooling
171            qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
172
173            ! --- Energy Budget of the leads (qlead, J.m-2) --- !
174            !     qlead is the energy received from the atm. in the leads.
175            !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2)
176            !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2)
177            IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN
178               ! upper bound for fhld: fhld should be equal to zqld
179               !                        but we have to make sure that this heat will not make the sst drop below the freezing point
180               !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos
181               !                        The following formulation is ok for both normal conditions and supercooling
182               fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
183                  &                                 - qsb_ice_bot(ji,jj) )
184               qlead(ji,jj) = 0._wp
185            ELSE
186               fhld (ji,jj) = 0._wp
187               ! upper bound for qlead: qlead should be equal to zqld
188               !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point.
189               !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb)
190               !                        and freezing point is reached if zqfr = zqld - qsb*a/dt
191               !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr
192               !                        The following formulation is ok for both normal conditions and supercooling
193               qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
194            ENDIF
195            !
196            ! If ice is landfast and ice concentration reaches its max
197            ! => stop ice formation in open water
198            IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp
199            !
200            ! If the grid cell is almost fully covered by ice (no leads)
201            ! => stop ice formation in open water
202            IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp
203            !
204            ! If ln_leadhfx is false
205            ! => do not use energy of the leads to melt sea-ice
206            IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp
207            !
208         END DO
209      END DO
210     
211      ! In case we bypass open-water ice formation
212      IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp
213      ! In case we bypass growing/melting from top and bottom
214      IF( .NOT. ln_icedH ) THEN
215         qsb_ice_bot(:,:) = 0._wp
216         fhld       (:,:) = 0._wp
217      ENDIF
218
219      !-------------------------------------------------------------------------------------------!
220      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
221      !-------------------------------------------------------------------------------------------!
222      DO jl = 1, jpl
223
224         ! select ice covered grid points
225         npti = 0 ; nptidx(:) = 0
226         DO jj = 1, jpj
227            DO ji = 1, jpi
228               IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
229                  npti         = npti  + 1
230                  nptidx(npti) = (jj - 1) * jpi + ji
231               ENDIF
232            END DO
233         END DO
234
235         IF( npti > 0 ) THEN  ! If there is no ice, do nothing.
236            !                                                               
237                              CALL ice_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- !
238            !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
239            !
240            s_i_new   (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp   ! --- some init --- !  (important to have them here)
241            dh_i_sum  (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm  (1:npti) = 0._wp 
242            dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp
243            dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
244            !                                     
245                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- !
246            !
247            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- !
248                              CALL ice_thd_dh                           ! Ice-Snow thickness   
249                              CALL ice_thd_pnd                          ! Melt ponds formation
250                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping
251            ENDIF
252                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !   
253            !
254                              CALL ice_thd_temp                     ! --- Temperature update --- !
255            !
256            IF( ln_icedH .AND. ln_virtual_itd ) &
257               &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- !
258            !
259            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- !
260            !
261                              CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
262            !                                                       ! --- & Move to 2D arrays --- !
263         ENDIF
264         !
265      END DO
266      !
267      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
268      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
269      !                   
270      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
271      !
272      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- !
273      !
274                              CALL ice_cor( kt , 2 )                ! --- Corrections --- !
275      !
276      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation     
277      !
278      ! convergence tests
279      IF( ln_zdf_chkcvg ) THEN
280         CALL iom_put( 'tice_cvgerr', ztice_cvgerr ) ; DEALLOCATE( ztice_cvgerr )
281         CALL iom_put( 'tice_cvgstp', ztice_cvgstp ) ; DEALLOCATE( ztice_cvgstp )
282      ENDIF
283      !
284      ! controls
285      IF( ln_icectl )   CALL ice_prt    (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints
286      IF( ln_ctl    )   CALL ice_prt3D  ('icethd')                                        ! prints
287      IF( ln_timing )   CALL timing_stop('icethd')                                        ! timing
288      !
289   END SUBROUTINE ice_thd 
290
291 
292   SUBROUTINE ice_thd_temp
293      !!-----------------------------------------------------------------------
294      !!                   ***  ROUTINE ice_thd_temp ***
295      !!                 
296      !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy
297      !!
298      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
299      !!-------------------------------------------------------------------
300      INTEGER  ::   ji, jk   ! dummy loop indices
301      REAL(wp) ::   ztmelts, zbbb, zccc  ! local scalar
302      !!-------------------------------------------------------------------
303      ! Recover ice temperature
304      DO jk = 1, nlay_i
305         DO ji = 1, npti
306            ztmelts       = -rTmlt * sz_i_1d(ji,jk)
307            ! Conversion q(S,T) -> T (second order equation)
308            zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus
309            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) )
310            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi
311           
312            ! mask temperature
313            rswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
314            t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
315         END DO
316      END DO 
317      !
318   END SUBROUTINE ice_thd_temp
319
320
321   SUBROUTINE ice_thd_mono
322      !!-----------------------------------------------------------------------
323      !!                   ***  ROUTINE ice_thd_mono ***
324      !!                 
325      !! ** Purpose :   Lateral melting in case virtual_itd
326      !!                          ( dA = A/2h dh )
327      !!-----------------------------------------------------------------------
328      INTEGER  ::   ji                 ! dummy loop indices
329      REAL(wp) ::   zhi_bef            ! ice thickness before thermo
330      REAL(wp) ::   zdh_mel, zda_mel   ! net melting
331      REAL(wp) ::   zvi, zvs           ! ice/snow volumes
332      !!-----------------------------------------------------------------------
333      !
334      DO ji = 1, npti
335         zdh_mel = MIN( 0._wp, dh_i_itm(ji) + dh_i_sum(ji) + dh_i_bom(ji) + dh_snowice(ji) + dh_i_sub(ji) )
336         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN
337            zvi          = a_i_1d(ji) * h_i_1d(ji)
338            zvs          = a_i_1d(ji) * h_s_1d(ji)
339            ! lateral melting = concentration change
340            zhi_bef     = h_i_1d(ji) - zdh_mel
341            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
342            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
343            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
344            ! adjust thickness
345            h_i_1d(ji) = zvi / a_i_1d(ji)           
346            h_s_1d(ji) = zvs / a_i_1d(ji)           
347            ! retrieve total concentration
348            at_i_1d(ji) = a_i_1d(ji)
349         END IF
350      END DO
351      !
352   END SUBROUTINE ice_thd_mono
353
354
355   SUBROUTINE ice_thd_1d2d( kl, kn )
356      !!-----------------------------------------------------------------------
357      !!                   ***  ROUTINE ice_thd_1d2d ***
358      !!                 
359      !! ** Purpose :   move arrays from 1d to 2d and the reverse
360      !!-----------------------------------------------------------------------
361      INTEGER, INTENT(in) ::   kl   ! index of the ice category
362      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
363      !
364      INTEGER ::   jk   ! dummy loop indices
365      !!-----------------------------------------------------------------------
366      !
367      SELECT CASE( kn )
368      !                    !---------------------!
369      CASE( 1 )            !==  from 2D to 1D  ==!
370         !                 !---------------------!
371         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
372         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
373         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
374         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
375         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
376         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
377         DO jk = 1, nlay_s
378            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
379            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
380         END DO
381         DO jk = 1, nlay_i
382            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
383            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
384            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
385         END DO
386         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
387         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
388         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) )
389         !
390         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            )
391         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d    (1:npti), qsr_ice (:,:,kl)     )
392         CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice (:,:,kl)     )
393         CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d   (1:npti), evap_ice(:,:,kl)     )
394         CALL tab_2d_1d( npti, nptidx(1:npti), dqns_ice_1d   (1:npti), dqns_ice(:,:,kl)     )
395         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d       (1:npti), t_bo                 )
396         CALL tab_2d_1d( npti, nptidx(1:npti), sprecip_1d    (1:npti), sprecip              ) 
397         CALL tab_2d_1d( npti, nptidx(1:npti), qsb_ice_bot_1d(1:npti), qsb_ice_bot          )
398         CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d       (1:npti), fhld                 )
399         
400         CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d    (1:npti), qml_ice    (:,:,kl) )
401         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d    (1:npti), qcn_ice    (:,:,kl) )
402         CALL tab_2d_1d( npti, nptidx(1:npti), qtr_ice_top_1d(1:npti), qtr_ice_top(:,:,kl) )
403         !
404         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni   )
405         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum   )
406         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sub_1d    (1:npti), wfx_sub       )
407         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub   )
408         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub   )
409         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub   )
410         !
411         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog          )
412         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom          )
413         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum          )
414         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni          )
415         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res          )
416         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          )
417         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          )
418         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          )
419         !
420         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          )
421         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom          )
422         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum          )
423         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni          )
424         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri          )
425         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res          )
426         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub          )
427         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam          )
428         !
429         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd       )
430         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_spr_1d    (1:npti), hfx_spr       )
431         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sum_1d    (1:npti), hfx_sum       )
432         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bom_1d    (1:npti), hfx_bom       )
433         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bog_1d    (1:npti), hfx_bog       )
434         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_lam_1d    (1:npti), hfx_lam       )
435         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dif_1d    (1:npti), hfx_dif       )
436         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d    (1:npti), hfx_opw       )
437         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_snw_1d    (1:npti), hfx_snw       )
438         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sub_1d    (1:npti), hfx_sub       )
439         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       )
440         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   )
441         !
442         ! ocean surface fields
443         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
444         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
445         CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m )
446         !
447         ! to update ice age
448         CALL tab_2d_1d( npti, nptidx(1:npti), o_i_1d (1:npti), o_i (:,:,kl) )
449         CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
450         !
451         ! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
452         DO jk = 1, nlay_i
453            WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
454         END DO
455         DO jk = 1, nlay_s
456            WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
457         END DO
458         !
459         !                 !---------------------!
460      CASE( 2 )            !==  from 1D to 2D  ==!
461         !                 !---------------------!
462         ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
463         DO jk = 1, nlay_i
464            e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * h_i_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_i
465         END DO
466         DO jk = 1, nlay_s
467            e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) * h_s_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_s
468         END DO
469         !
470         ! Change thickness to volume (replaces routine ice_var_eqv2glo)
471         v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti)
472         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti)
473         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
474         v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
475         v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
476         oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti)
477         
478         CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
479         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
480         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
481         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
482         CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
483         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
484         DO jk = 1, nlay_s
485            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
486            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
487         END DO
488         DO jk = 1, nlay_i
489            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
490            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
491            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
492         END DO
493         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
494         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
495         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) )
496         !
497         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
498         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
499         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sub_1d    (1:npti), wfx_sub     )
500         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub )
501         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub )
502         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub )
503         !
504         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog        )
505         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom        )
506         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        )
507         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni        )
508         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res        )
509         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        )
510         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        )
511         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
512         !
513         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        )
514         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom        )
515         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum        )
516         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni        )
517         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri        )
518         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res        )
519         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub        )
520         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam        )
521         !
522         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd     )
523         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d    (1:npti), hfx_spr     )
524         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d    (1:npti), hfx_sum     )
525         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d    (1:npti), hfx_bom     )
526         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d    (1:npti), hfx_bog     )
527         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_lam_1d    (1:npti), hfx_lam     )
528         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d    (1:npti), hfx_dif     )
529         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d    (1:npti), hfx_opw     )
530         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d    (1:npti), hfx_snw     )
531         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d    (1:npti), hfx_sub     )
532         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     )
533         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
534         !
535         CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) )
536         CALL tab_1d_2d( npti, nptidx(1:npti), qtr_ice_bot_1d(1:npti), qtr_ice_bot(:,:,kl) )
537         ! effective conductivity and 1st layer temperature (ln_cndflx=T)
538         CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) )
539         CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) )
540         ! SIMIP diagnostics         
541         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) )
542         CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot(:,:,kl) )
543         CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top(:,:,kl) )
544         ! extensive variables
545         CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
546         CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) )
547         CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
548         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
549         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) )
550         CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
551         ! check convergence of heat diffusion scheme
552         IF( ln_zdf_chkcvg ) THEN
553            CALL tab_1d_2d( npti, nptidx(1:npti), tice_cvgerr_1d(1:npti), ztice_cvgerr(:,:,kl) )
554            CALL tab_1d_2d( npti, nptidx(1:npti), tice_cvgstp_1d(1:npti), ztice_cvgstp(:,:,kl) )
555         ENDIF
556         !
557      END SELECT
558      !
559   END SUBROUTINE ice_thd_1d2d
560
561
562   SUBROUTINE ice_thd_init
563      !!-------------------------------------------------------------------
564      !!                   ***  ROUTINE ice_thd_init ***
565      !!                 
566      !! ** Purpose :   Physical constants and parameters associated with
567      !!                ice thermodynamics
568      !!
569      !! ** Method  :   Read the namthd namelist and check the parameters
570      !!                called at the first timestep (nit000)
571      !!
572      !! ** input   :   Namelist namthd
573      !!-------------------------------------------------------------------
574      INTEGER  ::   ios   ! Local integer output status for namelist read
575      !!
576      NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS, ln_leadhfx
577      !!-------------------------------------------------------------------
578      !
579      REWIND( numnam_ice_ref )              ! Namelist namthd in reference namelist : Ice thermodynamics
580      READ  ( numnam_ice_ref, namthd, IOSTAT = ios, ERR = 901)
581901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd in reference namelist' )
582      REWIND( numnam_ice_cfg )              ! Namelist namthd in configuration namelist : Ice thermodynamics
583      READ  ( numnam_ice_cfg, namthd, IOSTAT = ios, ERR = 902 )
584902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd in configuration namelist' )
585      IF(lwm) WRITE( numoni, namthd )
586      !
587      IF(lwp) THEN                          ! control print
588         WRITE(numout,*)
589         WRITE(numout,*) 'ice_thd_init: Ice Thermodynamics'
590         WRITE(numout,*) '~~~~~~~~~~~~'
591         WRITE(numout,*) '   Namelist namthd:'
592         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)                ln_icedH   = ', ln_icedH
593         WRITE(numout,*) '      activate lateral melting (T) or not (F)                              ln_icedA   = ', ln_icedA
594         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)                     ln_icedO   = ', ln_icedO
595         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)                ln_icedS   = ', ln_icedS
596         WRITE(numout,*) '      heat in the leads is used to melt sea-ice before warming the ocean   ln_leadhfx = ', ln_leadhfx
597     ENDIF
598      !
599                       CALL ice_thd_zdf_init   ! set ice heat diffusion parameters
600      IF( ln_icedA )   CALL ice_thd_da_init    ! set ice lateral melting parameters
601      IF( ln_icedO )   CALL ice_thd_do_init    ! set ice growth in open water parameters
602                       CALL ice_thd_sal_init   ! set ice salinity parameters
603                       CALL ice_thd_pnd_init   ! set melt ponds parameters
604      !
605   END SUBROUTINE ice_thd_init
606
607#else
608   !!----------------------------------------------------------------------
609   !!   Default option         Dummy module          NO  SI3 sea-ice model
610   !!----------------------------------------------------------------------
611#endif
612
613   !!======================================================================
614END MODULE icethd
Note: See TracBrowser for help on using the repository browser.