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

source: NEMO/trunk/src/ICE/iceistate.F90 @ 14088

Last change on this file since 14088 was 14086, checked in by cetlod, 4 years ago

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