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

Last change on this file was 15530, checked in by clem, 2 years ago

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