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/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE/icethd.F90 @ 11954

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

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.

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