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/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icethd.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 33.4 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, qns_tot, qsr_tot, 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 icetab         ! sea-ice: 1D <==> 2D transformation
33   USE icevar         ! sea-ice: operations
34   USE icectl         ! sea-ice: control print
35   !
36   USE in_out_manager ! I/O manager
37   USE lib_mpp        ! MPP library
38   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
39   USE lbclnk         ! lateral boundary conditions (or mpp links)
40   USE timing         ! Timing
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   ice_thd         ! called by limstp module
46   PUBLIC   ice_thd_init    ! called by ice_init
47
48   !!** namelist (namthd) **
49   LOGICAL ::   ln_icedH         ! activate ice thickness change from growing/melting (T) or not (F)
50   LOGICAL ::   ln_icedA         ! activate lateral melting param. (T) or not (F)
51   LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F)
52   LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
53
54   !! * Substitutions
55#  include "do_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
58   !! $Id$
59   !! Software governed by the CeCILL license (see ./LICENSE)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE ice_thd( kt )
64      !!-------------------------------------------------------------------
65      !!                ***  ROUTINE ice_thd  ***       
66      !! 
67      !! ** Purpose : This routine manages ice thermodynamics
68      !!         
69      !! ** Action : - computation of oceanic sensible heat flux at the ice base
70      !!                              energy budget in the leads
71      !!                              net fluxes on top of ice and of ocean
72      !!             - selection of grid cells with ice
73      !!                - call ice_thd_zdf  for vertical heat diffusion
74      !!                - call ice_thd_dh   for vertical ice growth and melt
75      !!                - call ice_thd_pnd  for melt ponds
76      !!                - call ice_thd_ent  for enthalpy remapping
77      !!                - call ice_thd_sal  for ice desalination
78      !!                - call ice_thd_temp to  retrieve temperature from ice enthalpy
79      !!                - call ice_thd_mono for extra lateral ice melt if active virtual thickness distribution
80      !!                - call ice_thd_da   for lateral ice melt
81      !!             - back to the geographic grid
82      !!                - call ice_thd_rem  for remapping thickness distribution
83      !!                - call ice_thd_do   for ice growth in leads
84      !!-------------------------------------------------------------------
85      INTEGER, INTENT(in) :: kt    ! number of iteration
86      !
87      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
88      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
89      REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04)
90      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient
91      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
92      !
93      !!-------------------------------------------------------------------
94      ! controls
95      IF( ln_timing    )   CALL timing_start('icethd')                                                             ! timing
96      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
97      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
98
99      IF( kt == nit000 .AND. lwp ) THEN
100         WRITE(numout,*)
101         WRITE(numout,*) 'ice_thd: sea-ice thermodynamics'
102         WRITE(numout,*) '~~~~~~~'
103      ENDIF
104     
105      !---------------------------------------------!
106      ! computation of friction velocity at T points
107      !---------------------------------------------!
108      IF( ln_icedyn ) THEN
109         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
110         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
111         DO_2D_00_00
112            zfric(ji,jj) = rn_cio * ( 0.5_wp *  &
113               &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
114               &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
115         END_2D
116      ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean
117         DO_2D_00_00
118            zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp *  &
119               &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
120               &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
121         END_2D
122      ENDIF
123      CALL lbc_lnk( 'icethd', zfric, 'T',  1. )
124      !
125      !--------------------------------------------------------------------!
126      ! Partial computation of forcing for the thermodynamic sea ice model
127      !--------------------------------------------------------------------!
128      DO_2D_11_11
129         rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
130         !
131         !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
132         !           !  practically no "direct lateral ablation"
133         !           
134         !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean
135         !           !  temperature and turbulent mixing (McPhee, 1992)
136         !
137         ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- !
138         zqld =  tmask(ji,jj,1) * rdt_ice *  &
139            &    ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) +  &
140            &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
141
142         ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- !
143         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)
144         zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0
145
146         ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2)
147         zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
148         qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
149
150         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
151         ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach
152         !                              the freezing point, so that we do not have SST < T_freeze
153         !                              This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
154
155         !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice
156         qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
157
158         ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting
159         ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget
160         IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN
161            fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
162            qlead(ji,jj) = 0._wp
163         ELSE
164            fhld (ji,jj) = 0._wp
165         ENDIF
166         !
167         ! Net heat flux on top of the ice-ocean [W.m-2]
168         ! ---------------------------------------------
169         qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 
170      END_2D
171     
172      ! In case we bypass open-water ice formation
173      IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp
174      ! In case we bypass growing/melting from top and bottom
175      IF( .NOT. ln_icedH ) THEN
176         qsb_ice_bot(:,:) = 0._wp
177         fhld       (:,:) = 0._wp
178      ENDIF
179
180      ! ---------------------------------------------------------------------
181      ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
182      ! ---------------------------------------------------------------------
183      !     First  step here              :  non solar + precip - qlead - qsensible
184      !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema)
185      !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar
186      qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean               
187         &             - qlead(:,:) * r1_rdtice                                &  ! heat flux taken from the ocean where there is open water ice formation
188         &             - at_i (:,:) * qsb_ice_bot(:,:)                         &  ! heat flux taken by sensible flux
189         &             - at_i (:,:) * fhld       (:,:)                            ! heat flux taken during bottom growth/melt
190      !                                                                           !    (fhld should be 0 while bott growth)
191      !-------------------------------------------------------------------------------------------!
192      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
193      !-------------------------------------------------------------------------------------------!
194      DO jl = 1, jpl
195
196         ! select ice covered grid points
197         npti = 0 ; nptidx(:) = 0
198         DO_2D_11_11
199            IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
200               npti         = npti  + 1
201               nptidx(npti) = (jj - 1) * jpi + ji
202            ENDIF
203         END_2D
204
205         IF( npti > 0 ) THEN  ! If there is no ice, do nothing.
206            !                                                               
207                              CALL ice_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- !
208            !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
209            !
210            s_i_new   (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp  ! --- some init --- !  (important to have them here)
211            dh_i_sum  (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm  (1:npti) = 0._wp 
212            dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp
213            dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
214            !                                     
215                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- !
216            !
217            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- !
218                              CALL ice_thd_dh                           ! Ice-Snow thickness   
219                              CALL ice_thd_pnd                          ! Melt ponds formation
220                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping
221            ENDIF
222                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !   
223            !
224                              CALL ice_thd_temp                     ! --- Temperature update --- !
225            !
226            IF( ln_icedH .AND. ln_virtual_itd ) &
227               &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- !
228            !
229            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- !
230            !
231                              CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
232            !                                                       ! --- & Move to 2D arrays --- !
233         ENDIF
234         !
235      END DO
236      !
237      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
238      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
239      !                   
240      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
241      !
242      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- !
243      !
244      ! controls
245      IF( ln_icectl )   CALL ice_prt    (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints
246      IF( sn_cfctl%l_prtctl )   &
247        &               CALL ice_prt3D  ('icethd')                                        ! prints
248      IF( ln_timing )   CALL timing_stop('icethd')                                        ! timing
249      !
250   END SUBROUTINE ice_thd 
251
252 
253   SUBROUTINE ice_thd_temp
254      !!-----------------------------------------------------------------------
255      !!                   ***  ROUTINE ice_thd_temp ***
256      !!                 
257      !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy
258      !!
259      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
260      !!-------------------------------------------------------------------
261      INTEGER  ::   ji, jk   ! dummy loop indices
262      REAL(wp) ::   ztmelts, zbbb, zccc  ! local scalar
263      !!-------------------------------------------------------------------
264      ! Recover ice temperature
265      DO jk = 1, nlay_i
266         DO ji = 1, npti
267            ztmelts       = -rTmlt * sz_i_1d(ji,jk)
268            ! Conversion q(S,T) -> T (second order equation)
269            zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus
270            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) )
271            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi
272           
273            ! mask temperature
274            rswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
275            t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
276         END DO
277      END DO 
278      !
279   END SUBROUTINE ice_thd_temp
280
281
282   SUBROUTINE ice_thd_mono
283      !!-----------------------------------------------------------------------
284      !!                   ***  ROUTINE ice_thd_mono ***
285      !!                 
286      !! ** Purpose :   Lateral melting in case virtual_itd
287      !!                          ( dA = A/2h dh )
288      !!-----------------------------------------------------------------------
289      INTEGER  ::   ji                 ! dummy loop indices
290      REAL(wp) ::   zhi_bef            ! ice thickness before thermo
291      REAL(wp) ::   zdh_mel, zda_mel   ! net melting
292      REAL(wp) ::   zvi, zvs           ! ice/snow volumes
293      !!-----------------------------------------------------------------------
294      !
295      DO ji = 1, npti
296         zdh_mel = MIN( 0._wp, dh_i_itm(ji) + dh_i_sum(ji) + dh_i_bom(ji) + dh_snowice(ji) + dh_i_sub(ji) )
297         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN
298            zvi          = a_i_1d(ji) * h_i_1d(ji)
299            zvs          = a_i_1d(ji) * h_s_1d(ji)
300            ! lateral melting = concentration change
301            zhi_bef     = h_i_1d(ji) - zdh_mel
302            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
303            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
304            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
305            ! adjust thickness
306            h_i_1d(ji) = zvi / a_i_1d(ji)           
307            h_s_1d(ji) = zvs / a_i_1d(ji)           
308            ! retrieve total concentration
309            at_i_1d(ji) = a_i_1d(ji)
310         END IF
311      END DO
312      !
313   END SUBROUTINE ice_thd_mono
314
315
316   SUBROUTINE ice_thd_1d2d( kl, kn )
317      !!-----------------------------------------------------------------------
318      !!                   ***  ROUTINE ice_thd_1d2d ***
319      !!                 
320      !! ** Purpose :   move arrays from 1d to 2d and the reverse
321      !!-----------------------------------------------------------------------
322      INTEGER, INTENT(in) ::   kl   ! index of the ice category
323      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
324      !
325      INTEGER ::   jk   ! dummy loop indices
326      !!-----------------------------------------------------------------------
327      !
328      SELECT CASE( kn )
329      !                    !---------------------!
330      CASE( 1 )            !==  from 2D to 1D  ==!
331         !                 !---------------------!
332         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
333         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
334         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
335         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
336         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
337         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
338         DO jk = 1, nlay_s
339            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
340            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
341         END DO
342         DO jk = 1, nlay_i
343            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
344            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
345            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
346         END DO
347         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
348         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
349         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
350         !
351         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            )
352         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d    (1:npti), qsr_ice (:,:,kl)     )
353         CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice (:,:,kl)     )
354         CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d   (1:npti), evap_ice(:,:,kl)     )
355         CALL tab_2d_1d( npti, nptidx(1:npti), dqns_ice_1d   (1:npti), dqns_ice(:,:,kl)     )
356         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d       (1:npti), t_bo                 )
357         CALL tab_2d_1d( npti, nptidx(1:npti), sprecip_1d    (1:npti), sprecip              ) 
358         CALL tab_2d_1d( npti, nptidx(1:npti), qsb_ice_bot_1d(1:npti), qsb_ice_bot          )
359         CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d       (1:npti), fhld                 )
360         
361         CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d    (1:npti), qml_ice    (:,:,kl) )
362         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d    (1:npti), qcn_ice    (:,:,kl) )
363         CALL tab_2d_1d( npti, nptidx(1:npti), qtr_ice_top_1d(1:npti), qtr_ice_top(:,:,kl) )
364         !
365         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni   )
366         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum   )
367         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sub_1d    (1:npti), wfx_sub       )
368         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub   )
369         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub   )
370         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub   )
371         !
372         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog          )
373         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom          )
374         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum          )
375         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni          )
376         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res          )
377         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          )
378         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          )
379         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          )
380         !
381         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          )
382         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom          )
383         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum          )
384         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni          )
385         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri          )
386         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res          )
387         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub          )
388         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam          )
389         !
390         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd       )
391         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_spr_1d    (1:npti), hfx_spr       )
392         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sum_1d    (1:npti), hfx_sum       )
393         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bom_1d    (1:npti), hfx_bom       )
394         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bog_1d    (1:npti), hfx_bog       )
395         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dif_1d    (1:npti), hfx_dif       )
396         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d    (1:npti), hfx_opw       )
397         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_snw_1d    (1:npti), hfx_snw       )
398         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sub_1d    (1:npti), hfx_sub       )
399         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       )
400         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   )
401         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem   )
402         CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     )
403         !
404         ! ocean surface fields
405         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
406         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
407         !
408         ! to update ice age
409         CALL tab_2d_1d( npti, nptidx(1:npti), o_i_1d (1:npti), o_i (:,:,kl) )
410         CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
411         !
412         ! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
413         DO jk = 1, nlay_i
414            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
415         END DO
416         DO jk = 1, nlay_s
417            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
418         END DO
419         !
420         !                 !---------------------!
421      CASE( 2 )            !==  from 1D to 2D  ==!
422         !                 !---------------------!
423         ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
424         DO jk = 1, nlay_i
425            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
426         END DO
427         DO jk = 1, nlay_s
428            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
429         END DO
430         !
431         ! Change thickness to volume (replaces routine ice_var_eqv2glo)
432         v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti)
433         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti)
434         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
435         v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
436         oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti)
437         
438         CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
439         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
440         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
441         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
442         CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
443         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
444         DO jk = 1, nlay_s
445            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
446            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
447         END DO
448         DO jk = 1, nlay_i
449            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
450            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
451            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
452         END DO
453         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
454         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
455         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
456         !
457         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
458         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
459         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sub_1d    (1:npti), wfx_sub     )
460         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub )
461         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub )
462         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub )
463         !
464         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog        )
465         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom        )
466         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        )
467         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni        )
468         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res        )
469         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        )
470         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        )
471         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
472         !
473         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        )
474         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom        )
475         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum        )
476         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni        )
477         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri        )
478         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res        )
479         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub        )
480         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam        )
481         !
482         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd     )
483         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d    (1:npti), hfx_spr     )
484         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d    (1:npti), hfx_sum     )
485         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d    (1:npti), hfx_bom     )
486         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d    (1:npti), hfx_bog     )
487         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d    (1:npti), hfx_dif     )
488         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d    (1:npti), hfx_opw     )
489         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d    (1:npti), hfx_snw     )
490         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d    (1:npti), hfx_sub     )
491         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     )
492         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
493         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )
494         CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   )
495         !
496         CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) )
497         CALL tab_1d_2d( npti, nptidx(1:npti), qtr_ice_bot_1d(1:npti), qtr_ice_bot(:,:,kl) )
498         ! effective conductivity and 1st layer temperature (ln_cndflx=T)
499         CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) )
500         CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) )
501         ! SIMIP diagnostics         
502         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) )
503         CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot(:,:,kl) )
504         CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top(:,:,kl) )
505         ! extensive variables
506         CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
507         CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) )
508         CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
509         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
510         CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
511         !
512      END SELECT
513      !
514   END SUBROUTINE ice_thd_1d2d
515
516
517   SUBROUTINE ice_thd_init
518      !!-------------------------------------------------------------------
519      !!                   ***  ROUTINE ice_thd_init ***
520      !!                 
521      !! ** Purpose :   Physical constants and parameters associated with
522      !!                ice thermodynamics
523      !!
524      !! ** Method  :   Read the namthd namelist and check the parameters
525      !!                called at the first timestep (nit000)
526      !!
527      !! ** input   :   Namelist namthd
528      !!-------------------------------------------------------------------
529      INTEGER  ::   ios   ! Local integer output status for namelist read
530      !!
531      NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS
532      !!-------------------------------------------------------------------
533      !
534      READ  ( numnam_ice_ref, namthd, IOSTAT = ios, ERR = 901)
535901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd in reference namelist' )
536      READ  ( numnam_ice_cfg, namthd, IOSTAT = ios, ERR = 902 )
537902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd in configuration namelist' )
538      IF(lwm) WRITE( numoni, namthd )
539      !
540      IF(lwp) THEN                          ! control print
541         WRITE(numout,*)
542         WRITE(numout,*) 'ice_thd_init: Ice Thermodynamics'
543         WRITE(numout,*) '~~~~~~~~~~~~'
544         WRITE(numout,*) '   Namelist namthd:'
545         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)   ln_icedH  = ', ln_icedH
546         WRITE(numout,*) '      activate lateral melting (T) or not (F)                 ln_icedA  = ', ln_icedA
547         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)        ln_icedO  = ', ln_icedO
548         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)   ln_icedS  = ', ln_icedS
549     ENDIF
550      !
551                       CALL ice_thd_zdf_init   ! set ice heat diffusion parameters
552      IF( ln_icedA )   CALL ice_thd_da_init    ! set ice lateral melting parameters
553      IF( ln_icedO )   CALL ice_thd_do_init    ! set ice growth in open water parameters
554                       CALL ice_thd_sal_init   ! set ice salinity parameters
555                       CALL ice_thd_pnd_init   ! set melt ponds parameters
556      !
557   END SUBROUTINE ice_thd_init
558
559#else
560   !!----------------------------------------------------------------------
561   !!   Default option         Dummy module          NO  SI3 sea-ice model
562   !!----------------------------------------------------------------------
563#endif
564
565   !!======================================================================
566END MODULE icethd
Note: See TracBrowser for help on using the repository browser.