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.
iceistate.F90 in NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/ICE/iceistate.F90 @ 14018

Last change on this file since 14018 was 13998, checked in by techene, 4 years ago

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