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 @ 14558

Last change on this file since 14558 was 14143, checked in by techene, 3 years ago

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