source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE/iceistate.F90 @ 14021

Last change on this file since 14021 was 14021, checked in by laurent, 5 months ago

Caught up with trunk rev 14020…

  • Property svn:keywords set to Id
File size: 30.2 KB
Line 
1MODULE iceistate
2   !!======================================================================
3   !!                     ***  MODULE  iceistate  ***
4   !!   sea-ice : Initialization of ice variables
5   !!======================================================================
6   !! History :  2.0  !  2004-01  (C. Ethe, G. Madec) Original code
7   !!            3.0  !  2007     (M. Vancoppenolle)  Rewrite for ice cats
8   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_istate       :  initialization of diagnostics ice variables
15   !!   ice_istate_init  :  initialization of ice state and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst         ! physical constant
18   USE oce            ! dynamics and tracers variables
19   USE dom_oce        ! ocean domain
20   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd
21   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
22   USE eosbn2         ! equation of state
23   USE domvvl         ! Variable volume
24   USE ice            ! sea-ice: variables
25   USE ice1D          ! sea-ice: thermodynamics variables
26   USE icetab         ! sea-ice: 1D <==> 2D transformation
27   USE icevar         ! sea-ice: operations
28   !
29   USE in_out_manager ! I/O manager
30   USE iom            ! I/O manager library
31   USE lib_mpp        ! MPP library
32   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
33   USE fldread        ! read input fields
34
35# if defined key_agrif
36   USE agrif_oce
37   USE agrif_ice
38   USE agrif_ice_interp
39# endif
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   ice_istate        ! called by icestp.F90
45   PUBLIC   ice_istate_init   ! called by icestp.F90
46   !
47   !                             !! ** namelist (namini) **
48   LOGICAL, PUBLIC  ::   ln_iceini        !: Ice initialization or not
49   INTEGER, PUBLIC  ::   nn_iceini_file   !: Ice initialization:
50                                  !        0 = Initialise sea ice based on SSTs
51                                  !        1 = Initialise sea ice from single category netcdf file
52                                  !        2 = Initialise sea ice from multi category restart file
53   REAL(wp) ::   rn_thres_sst
54   REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
55   REAL(wp) ::   rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
56   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n, rn_hld_ini_n
57   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s, rn_hld_ini_s
58   !
59   !                              ! if nn_iceini_file = 1
60   INTEGER , PARAMETER ::   jpfldi = 10          ! maximum number of files to read
61   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m)
62   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m)
63   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction     (-)
64   INTEGER , PARAMETER ::   jp_smi = 4           ! index of ice salinity     (g/kg)
65   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temperature  (K)
66   INTEGER , PARAMETER ::   jp_tsu = 6           ! index of ice surface temp (K)
67   INTEGER , PARAMETER ::   jp_tms = 7           ! index of snw temperature  (K)
68   INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-)
69   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
70   INTEGER , PARAMETER ::   jp_hld = 10          ! index of pnd lid depth    (m)
71   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
72
73   !! * Substitutions
74#  include "do_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
77   !! $Id$
78   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
79   !!----------------------------------------------------------------------
80CONTAINS
81
82   SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa )
83      !!-------------------------------------------------------------------
84      !!                    ***  ROUTINE ice_istate  ***
85      !!
86      !! ** Purpose :   defined the sea-ice initial state
87      !!
88      !! ** Method  :   This routine will put some ice where ocean
89      !!                is at the freezing point, then fill in ice
90      !!                state variables using prescribed initial
91      !!                values in the namelist
92      !!
93      !! ** Steps   :   1) Set initial surface and basal temperatures
94      !!                2) Recompute or read sea ice state variables
95      !!                3) Fill in space-dependent arrays for state variables
96      !!                4) snow-ice mass computation
97      !!
98      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
99      !!              where there is no ice
100      !!--------------------------------------------------------------------
101      INTEGER, INTENT(in) :: kt            ! time step
102      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
103      !
104      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
105      REAL(wp) ::   ztmelts
106      INTEGER , DIMENSION(4)           ::   itest
107      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
108      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
109      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
110      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
111      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini, zhlid_ini            !data from namelist or nc file
112      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
113      !!
114      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d
115      !--------------------------------------------------------------------
116
117      IF(lwp) WRITE(numout,*)
118      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
119      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
120
121      !---------------------------
122      ! 1) 1st init. of the fields
123      !---------------------------
124      !
125      ! basal temperature (considered at freezing point)   [Kelvin]
126      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
127      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
128      !
129      ! surface temperature and conductivity
130      DO jl = 1, jpl
131         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface
132         cnd_ice(:,:,jl) = 0._wp               ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
133      END DO
134      !
135      ! ice and snw temperatures
136      DO jl = 1, jpl
137         DO jk = 1, nlay_i
138            t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
139         END DO
140         DO jk = 1, nlay_s
141            t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
142         END DO
143      END DO
144      !
145      ! specific temperatures for coupled runs
146      tn_ice (:,:,:) = t_i (:,:,1,:)
147      t1_ice (:,:,:) = t_i (:,:,1,:)
148
149      ! heat contents
150      e_i (:,:,:,:) = 0._wp
151      e_s (:,:,:,:) = 0._wp
152
153      ! general fields
154      a_i (:,:,:) = 0._wp
155      v_i (:,:,:) = 0._wp
156      v_s (:,:,:) = 0._wp
157      sv_i(:,:,:) = 0._wp
158      oa_i(:,:,:) = 0._wp
159      !
160      h_i (:,:,:) = 0._wp
161      h_s (:,:,:) = 0._wp
162      s_i (:,:,:) = 0._wp
163      o_i (:,:,:) = 0._wp
164      !
165      ! melt ponds
166      a_ip     (:,:,:) = 0._wp
167      v_ip     (:,:,:) = 0._wp
168      v_il     (:,:,:) = 0._wp
169      a_ip_eff (:,:,:) = 0._wp
170      h_ip     (:,:,:) = 0._wp
171      h_il     (:,:,:) = 0._wp
172      !
173      ! ice velocities
174      u_ice (:,:) = 0._wp
175      v_ice (:,:) = 0._wp
176      !
177      !------------------------------------------------------------------------
178      ! 2) overwrite some of the fields with namelist parameters or netcdf file
179      !------------------------------------------------------------------------
180      IF( ln_iceini ) THEN
181         !
182         IF( Agrif_Root() ) THEN
183            !                             !---------------!
184            IF( nn_iceini_file == 1 )THEN ! Read a file   !
185               !                          !---------------!
186               WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
187               ELSEWHERE                     ;   zswitch(:,:) = 0._wp
188               END WHERE
189               !
190               CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
191               !
192               ! -- mandatory fields -- !
193               zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
194               zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
195               zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
196
197               ! -- optional fields -- !
198               !    if fields do not exist then set them to the values present in the namelist (except for temperatures)
199               !
200               ! ice salinity
201               IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
202                  &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
203               !
204               ! temperatures
205               IF    ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
206                  &    TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
207                  si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
208                  si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
209                  si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
210               ENDIF
211               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
212                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
213               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
214                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
215               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s
216                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
217               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i
218                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
219               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su
220                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
221               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i
222                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
223               !
224               ! pond concentration
225               IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
226                  &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
227                  &                              * si(jp_ati)%fnow(:,:,1)
228               !
229               ! pond depth
230               IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
231                  &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
232               !
233               ! pond lid depth
234               IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) &
235                  &     si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
236               !
237               zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
238               ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
239               zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
240               ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
241               zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
242               zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
243               zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1)
244               !
245               ! change the switch for the following
246               WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1)
247               ELSEWHERE                         ;   zswitch(:,:) = 0._wp
248               END WHERE
249
250               !                          !---------------!
251            ELSE                          ! Read namelist !
252               !                          !---------------!
253               ! no ice if (sst - Tfreez) >= thresold
254               WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp
255               ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
256               END WHERE
257               !
258               ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
259               WHERE( ff_t(:,:) >= 0._wp )
260                  zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
261                  zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
262                  zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
263                  zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
264                  ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
265                  zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
266                  ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
267                  zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
268                  zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
269                  zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:)
270               ELSEWHERE
271                  zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
272                  zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
273                  zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
274                  zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
275                  ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
276                  zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
277                  ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
278                  zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
279                  zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
280                  zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:)
281               END WHERE
282               !
283            ENDIF
284
285
286
287            ! make sure ponds = 0 if no ponds scheme
288            IF ( .NOT.ln_pnd ) THEN
289               zapnd_ini(:,:) = 0._wp
290               zhpnd_ini(:,:) = 0._wp
291               zhlid_ini(:,:) = 0._wp
292            ENDIF
293
294            IF ( .NOT.ln_pnd_lids ) THEN
295               zhlid_ini(:,:) = 0._wp
296            ENDIF
297
298            !----------------!
299            ! 3) fill fields !
300            !----------------!
301            ! select ice covered grid points
302            npti = 0 ; nptidx(:) = 0
303            DO_2D( 1, 1, 1, 1 )
304               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
305                  npti         = npti  + 1
306                  nptidx(npti) = (jj - 1) * jpi + ji
307               ENDIF
308            END_2D
309
310            ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
311            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
312            CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
313            CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
314            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
315            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
316            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
317            CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
318            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
319            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
320            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_ini )
321
322            ! allocate temporary arrays
323            ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), &
324               &      zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), &
325               &      zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) )
326
327            ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
328            CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                  &
329               &              zhi_2d          , zhs_2d          , zai_2d         ,                  &
330               &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti),                  &
331               &              s_i_1d(1:npti)  , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), &
332               &              zti_2d          , zts_2d          , ztsu_2d        ,                  &
333               &              zsi_2d          , zaip_2d         , zhip_2d        , zhil_2d )
334
335            ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
336            DO jl = 1, jpl
337               zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
338               zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
339            END DO
340            CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
341            CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
342            CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
343            CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
344            CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
345            CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
346            CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
347            CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
348            CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
349            CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d  , h_il   )
350
351            ! deallocate temporary arrays
352            DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
353               &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d )
354
355            ! calculate extensive and intensive variables
356            CALL ice_var_salprof ! for sz_i
357            DO jl = 1, jpl
358               DO_2D( 1, 1, 1, 1 )
359                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
360                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
361                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
362               END_2D
363            END DO
364            !
365            DO jl = 1, jpl
366               DO_3D( 1, 1, 1, 1, 1, nlay_s )
367                  t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
368                  e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
369                     &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
370               END_3D
371            END DO
372            !
373            DO jl = 1, jpl
374               DO_3D( 1, 1, 1, 1, 1, nlay_i )
375                  t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
376                  ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
377                  e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
378                     &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
379                     &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
380                     &                       - rcp   * ( ztmelts - rt0 ) )
381               END_3D
382            END DO
383
384#if  defined key_agrif
385         ELSE
386
387            Agrif_SpecialValue    = -9999.
388            Agrif_UseSpecialValue = .TRUE.
389            CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice)
390            use_sign_north = .TRUE.
391            sign_north = -1.
392            CALL Agrif_init_variable(u_iceini_id  ,procname=interp_u_ice)
393            CALL Agrif_init_variable(v_iceini_id  ,procname=interp_v_ice)
394            Agrif_SpecialValue    = 0._wp
395            use_sign_north = .FALSE.
396            Agrif_UseSpecialValue = .FALSE.
397        ! lbc ????
398   ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, v_il, t_su, e_s, e_i
399            CALL ice_var_glo2eqv
400            CALL ice_var_zapsmall
401            CALL ice_var_agg(2)
402#endif
403         ENDIF ! Agrif_Root
404         !
405         ! Melt ponds
406         WHERE( a_i > epsi10 )   ;   a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
407         ELSEWHERE               ;   a_ip_eff(:,:,:) = 0._wp
408         END WHERE
409         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
410         v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
411
412         ! specific temperatures for coupled runs
413         tn_ice(:,:,:) = t_su(:,:,:)
414         t1_ice(:,:,:) = t_i (:,:,1,:)
415         !
416         ! ice concentration should not exceed amax
417         at_i(:,:) = SUM( a_i, dim=3 )
418         DO jl = 1, jpl
419            WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)
420         END DO
421         at_i(:,:) = SUM( a_i, dim=3 )
422         !
423      ENDIF ! ln_iceini
424      !
425      !----------------------------------------------
426      ! 4) Snow-ice mass (case ice is fully embedded)
427      !----------------------------------------------
428      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3  )   ! snow+ice mass
429      snwice_mass_b(:,:) = snwice_mass(:,:)
430      !
431      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
432         !
433         ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rho0
434         ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0
435         !
436         IF( .NOT.ln_linssh )   CALL dom_vvl_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
437! !!st
438!          IF( .NOT.ln_linssh ) THEN
439!             !
440!             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:)
441!             ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
442!             !
443!             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
444!                e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:)
445!                e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
446!                e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm)
447!             END DO
448!             !
449!             ! Reconstruction of all vertical scale factors at now and before time-steps
450!             ! =========================================================================
451!             ! Horizontal scale factor interpolations
452!             ! --------------------------------------
453!             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
454!             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
455!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
456!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
457!             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
458!             ! Vertical scale factor interpolations
459!             ! ------------------------------------
460!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
461!             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
462!             CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
463!             CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
464!             CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
465!             ! t- and w- points depth
466!             ! ----------------------
467!             !!gm not sure of that....
468!             gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
469!             gdepw(:,:,1,Kmm) = 0.0_wp
470!             gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
471!             DO jk = 2, jpk
472!                gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk  ,Kmm)
473!                gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
474!                gde3w(:,:,jk) = gdept(:,:,jk  ,Kmm) - ssh (:,:,Kmm)
475!             END DO
476!          ENDIF
477      ENDIF
478
479      !!clem: output of initial state should be written here but it is impossible because
480      !!      the ocean and ice are in the same file
481      !!      CALL dia_wri_state( 'output.init' )
482      !
483   END SUBROUTINE ice_istate
484
485
486   SUBROUTINE ice_istate_init
487      !!-------------------------------------------------------------------
488      !!                   ***  ROUTINE ice_istate_init  ***
489      !!
490      !! ** Purpose :   Definition of initial state of the ice
491      !!
492      !! ** Method  :   Read the namini namelist and check the parameter
493      !!              values called at the first timestep (nit000)
494      !!
495      !! ** input   :  Namelist namini
496      !!
497      !!-----------------------------------------------------------------------------
498      INTEGER ::   ios   ! Local integer output status for namelist read
499      INTEGER ::   ifpr, ierror
500      !
501      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
502      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld
503      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
504      !
505      NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, &
506         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
507         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
508         &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
509         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, &
510         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir
511      !!-----------------------------------------------------------------------------
512      !
513      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
514901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namini in reference namelist' )
515      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
516902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namini in configuration namelist' )
517      IF(lwm) WRITE ( numoni, namini )
518      !
519      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
520      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
521      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
522      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd   ;   slf_i(jp_hld) = sn_hld
523      !
524      IF(lwp) THEN                          ! control print
525         WRITE(numout,*)
526         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
527         WRITE(numout,*) '~~~~~~~~~~~~~~~'
528         WRITE(numout,*) '   Namelist namini:'
529         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini
530         WRITE(numout,*) '      ice initialization from a netcdf file            nn_iceini_file = ', nn_iceini_file
531         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
532         IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN
533            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s
534            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
535            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s
536            WRITE(numout,*) '      initial ice salinity  in the north-south         rn_smi_ini     = ', rn_smi_ini_n,rn_smi_ini_s
537            WRITE(numout,*) '      initial surf temperat in the north-south         rn_tsu_ini     = ', rn_tsu_ini_n,rn_tsu_ini_s
538            WRITE(numout,*) '      initial ice temperat  in the north-south         rn_tmi_ini     = ', rn_tmi_ini_n,rn_tmi_ini_s
539            WRITE(numout,*) '      initial snw temperat  in the north-south         rn_tms_ini     = ', rn_tms_ini_n,rn_tms_ini_s
540            WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
541            WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
542            WRITE(numout,*) '      initial pnd lid depth in the north-south         rn_hld_ini     = ', rn_hld_ini_n,rn_hld_ini_s
543         ENDIF
544      ENDIF
545      !
546      IF( nn_iceini_file == 1 ) THEN                      ! Ice initialization using input file
547         !
548         ! set si structure
549         ALLOCATE( si(jpfldi), STAT=ierror )
550         IF( ierror > 0 ) THEN
551            CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' )   ;   RETURN
552         ENDIF
553         !
554         DO ifpr = 1, jpfldi
555            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
556            IF( slf_i(ifpr)%ln_tint )  ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
557         END DO
558         !
559         ! fill si with slf_i and control print
560         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
561         !
562      ENDIF
563      !
564      IF( .NOT.ln_pnd ) THEN
565         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
566         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
567         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
568         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' )
569      ENDIF
570      !
571      IF( .NOT.ln_pnd_lids ) THEN
572         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
573      ENDIF
574      !
575   END SUBROUTINE ice_istate_init
576
577#else
578   !!----------------------------------------------------------------------
579   !!   Default option :         Empty module         NO SI3 sea-ice model
580   !!----------------------------------------------------------------------
581#endif
582
583   !!======================================================================
584END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.