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

source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceistate.F90 @ 13471

Last change on this file since 13471 was 13471, checked in by smasson, 4 years ago

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