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

source: NEMO/branches/2020/ticket2487/src/ICE/iceistate.F90 @ 13299

Last change on this file since 13299 was 13299, checked in by smueller, 3 years ago

Addition of a namelist parameter to enable the selectability of sea-level compensation for non-zero sea-ice/snow mass (see ticket #2487)
  • Property svn:keywords set to Id
File size: 28.7 KB
Line 
1MODULE iceistate
2   !!======================================================================
3   !!                     ***  MODULE  iceistate  ***
4   !!   sea-ice : Initialization of ice variables
5   !!======================================================================
6   !! History :  2.0  !  2004-01  (C. Ethe, G. Madec) Original code
7   !!            3.0  !  2007     (M. Vancoppenolle)  Rewrite for ice cats
8   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_istate       :  initialization of diagnostics ice variables
15   !!   ice_istate_init  :  initialization of ice state and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst         ! physical constant
18   USE oce            ! dynamics and tracers variables
19   USE dom_oce        ! ocean domain
20   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd , ln_ice_sladj
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   LOGICAL, PUBLIC  ::   ln_iceini_file   !: Ice initialization from 2D netcdf file
44   REAL(wp) ::   rn_thres_sst
45   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
46   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
47   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n
48   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s
49   !
50   !                              ! if ln_iceini_file = T
51   INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read
52   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m)
53   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m)
54   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction     (-)
55   INTEGER , PARAMETER ::   jp_smi = 4           ! index of ice salinity     (g/kg)
56   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temperature  (K)
57   INTEGER , PARAMETER ::   jp_tsu = 6           ! index of ice surface temp (K)
58   INTEGER , PARAMETER ::   jp_tms = 7           ! index of snw temperature  (K)
59   INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-)
60   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
62   !
63#if defined key_agrif
64   REAL(wp), PUBLIC ::   rsshadj          !: initial mean ssh adjustment due to initial ice+snow mass
65#endif
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 the ice thickness distribution using gaussian
88      !!                4) Fill in space-dependent arrays for state variables
89      !!                5) snow-ice mass computation
90      !!                6) store before fields
91      !!
92      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
93      !!              where there is no ice
94      !!--------------------------------------------------------------------
95      INTEGER, INTENT(in) ::   kt   ! time step
96      !!
97      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
98      REAL(wp) ::   ztmelts, z1_area, zsshadj
99      INTEGER , DIMENSION(4)           ::   itest
100      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
101      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
102      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
103      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
104      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file
105      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
106      !!
107      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
108      !--------------------------------------------------------------------
109
110      IF(lwp) WRITE(numout,*)
111      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
112      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
113
114      !---------------------------
115      ! 1) 1st init. of the fields
116      !---------------------------
117      !
118      ! basal temperature (considered at freezing point)   [Kelvin]
119      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
120      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
121      !
122      ! surface temperature and conductivity
123      DO jl = 1, jpl
124         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface
125         cnd_ice(:,:,jl) = 0._wp               ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
126      END DO
127      !
128      ! ice and snw temperatures
129      DO jl = 1, jpl
130         DO jk = 1, nlay_i
131            t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
132         END DO
133         DO jk = 1, nlay_s
134            t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
135         END DO
136      END DO
137      !
138      ! specific temperatures for coupled runs
139      tn_ice (:,:,:) = t_i (:,:,1,:)
140      t1_ice (:,:,:) = t_i (:,:,1,:)
141
142      ! heat contents
143      e_i (:,:,:,:) = 0._wp
144      e_s (:,:,:,:) = 0._wp
145     
146      ! general fields
147      a_i (:,:,:) = 0._wp
148      v_i (:,:,:) = 0._wp
149      v_s (:,:,:) = 0._wp
150      sv_i(:,:,:) = 0._wp
151      oa_i(:,:,:) = 0._wp
152      !
153      h_i (:,:,:) = 0._wp
154      h_s (:,:,:) = 0._wp
155      s_i (:,:,:) = 0._wp
156      o_i (:,:,:) = 0._wp
157      !
158      ! melt ponds
159      a_ip     (:,:,:) = 0._wp
160      v_ip     (:,:,:) = 0._wp
161      a_ip_frac(:,:,:) = 0._wp
162      h_ip     (:,:,:) = 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( ln_iceini_file )THEN      ! Read a file   !
174            !                          !---------------!
175            WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
176            ELSEWHERE                     ;   zswitch(:,:) = 0._wp
177            END WHERE
178            !
179            CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
180            !
181            ! -- mandatory fields -- !
182            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
183            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
184            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
185
186            ! -- optional fields -- !
187            !    if fields do not exist then set them to the values present in the namelist (except for temperatures)
188            !
189            ! ice salinity
190            IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
191               &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
192            !
193            ! temperatures
194            IF    ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
195               &    TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
196               si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
197               si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
198               si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
199            ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
200               si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
201            ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
202               si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
203            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s
204               si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
205            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i
206               si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
207            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su
208               si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
209            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i
210               si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
211            ENDIF
212            !
213            ! pond concentration
214            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
215               &     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.
216               &                              * si(jp_ati)%fnow(:,:,1) 
217            !
218            ! pond depth
219            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
220               &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
221            !
222            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
223            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
224            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
225            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
226            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
227            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
228            !
229            ! change the switch for the following
230            WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
231            ELSEWHERE                         ;   zswitch(:,:) = 0._wp
232            END WHERE
233            !                          !---------------!
234         ELSE                          ! Read namelist !
235            !                          !---------------!
236            ! no ice if (sst - Tfreez) >= thresold
237            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
238            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
239            END WHERE
240            !
241            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
242            WHERE( ff_t(:,:) >= 0._wp )
243               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
244               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
245               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
246               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
247               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
248               zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
249               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
250               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
251               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
252            ELSEWHERE
253               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
254               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
255               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
256               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
257               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
258               zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
259               ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
260               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
261               zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
262            END WHERE
263            !
264         ENDIF
265
266         ! make sure ponds = 0 if no ponds scheme
267         IF ( .NOT.ln_pnd ) THEN
268            zapnd_ini(:,:) = 0._wp
269            zhpnd_ini(:,:) = 0._wp
270         ENDIF
271         
272         !-------------!
273         ! fill fields !
274         !-------------!
275         ! select ice covered grid points
276         npti = 0 ; nptidx(:) = 0
277         DO jj = 1, jpj
278            DO ji = 1, jpi
279               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
280                  npti         = npti  + 1
281                  nptidx(npti) = (jj - 1) * jpi + ji
282               ENDIF
283            END DO
284         END DO
285
286         ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
287         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
288         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
289         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
290         CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
291         CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
292         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
293         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
294         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
295         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
296
297         ! allocate temporary arrays
298         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
299            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
300         
301         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
302         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
303            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   &
304            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), &
305            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d )
306
307         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
308         DO jl = 1, jpl
309            zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
310            zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
311         END DO
312         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
313         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
314         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
315         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
316         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
317         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
318         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
319         CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
320         CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
321
322         ! deallocate temporary arrays
323         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
324            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
325
326         ! calculate extensive and intensive variables
327         CALL ice_var_salprof ! for sz_i
328         DO jl = 1, jpl
329            DO jj = 1, jpj
330               DO ji = 1, jpi
331                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
332                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
333                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
334               END DO
335            END DO
336         END DO
337         !
338         DO jl = 1, jpl
339            DO jk = 1, nlay_s
340               DO jj = 1, jpj
341                  DO ji = 1, jpi
342                     t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
343                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
344                        &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
345                  END DO
346               END DO
347            END DO
348         END DO
349         !
350         DO jl = 1, jpl
351            DO jk = 1, nlay_i
352               DO jj = 1, jpj
353                  DO ji = 1, jpi
354                     t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
355                     ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
356                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
357                        &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
358                        &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
359                        &                       - rcp   * ( ztmelts - rt0 ) )
360                  END DO
361               END DO
362            END DO
363         END DO
364
365         ! Melt ponds
366         WHERE( a_i > epsi10 )
367            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
368         ELSEWHERE
369            a_ip_frac(:,:,:) = 0._wp
370         END WHERE
371         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
372         
373         ! specific temperatures for coupled runs
374         tn_ice(:,:,:) = t_su(:,:,:)
375         t1_ice(:,:,:) = t_i (:,:,1,:)
376         !
377      ENDIF ! ln_iceini
378      !
379      at_i(:,:) = SUM( a_i, dim=3 )
380      !
381      !-----------------
382      ! 3) Snow-ice mass
383      !-----------------
384      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass
385      snwice_mass_b(:,:) = snwice_mass(:,:)
386      !
387      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
388         !
389         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
390         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
391         !
392         IF( .NOT.ln_linssh ) THEN
393            !
394            z2d(:,:) = 1._wp + sshn(:,:) * ssmask(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) )
395            !
396            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
397               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
398               e3t_b(:,:,jk) = e3t_n(:,:,jk)
399               e3t_a(:,:,jk) = e3t_n(:,:,jk)
400            END DO
401            !
402            ! Reconstruction of all vertical scale factors at now and before time-steps
403            ! =========================================================================
404            ! Horizontal scale factor interpolations
405            ! --------------------------------------
406            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
407            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
408            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
409            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
410            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
411            ! Vertical scale factor interpolations
412            ! ------------------------------------
413            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
414            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
415            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
416            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
417            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
418            ! t- and w- points depth
419            ! ----------------------
420            !!gm not sure of that....
421            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
422            gdepw_n(:,:,1) = 0.0_wp
423            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
424            DO jk = 2, jpk
425               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
426               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
427               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
428            END DO
429         ENDIF
430      ELSEIF( ln_ice_sladj ) THEN       ! adjustment of sea level due to initial snow+ice mass
431         z1_area = 1.0_wp / glob_sum( 'iceistate', e1e2t(:,:) )
432         zsshadj = glob_sum( 'iceistate', e1e2t(:,:) * snwice_mass(:,:) ) * r1_rau0 * z1_area
433#if defined key_agrif
434         ! Override ssh adjustment in nested domains by the root-domain ssh
435         ! adjustment; store the adjustment value in a global module variable to
436         ! make it retrievable in nested domains
437         IF( .NOT. Agrif_Root() ) zsshadj = Agrif_Parent(rsshadj)
438         rsshadj = zsshadj
439#endif
440         IF(lwp) WRITE(numout,'(A35,F10.6,A21)') 'iceistate:   sea level adjusted by ', -1.0_wp * zsshadj, ' m to compensate for'
441         IF(lwp) WRITE(numout,*) '             the initial snow+ice mass'
442         sshn(:,:) = sshn(:,:) - zsshadj
443         sshb(:,:) = sshb(:,:) - zsshadj
444      ENDIF
445     
446      !------------------------------------
447      ! 4) store fields at before time-step
448      !------------------------------------
449      ! it is only necessary for the 1st interpolation by Agrif
450      a_i_b  (:,:,:)   = a_i  (:,:,:)
451      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
452      v_i_b  (:,:,:)   = v_i  (:,:,:)
453      v_s_b  (:,:,:)   = v_s  (:,:,:)
454      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
455      sv_i_b (:,:,:)   = sv_i (:,:,:)
456      oa_i_b (:,:,:)   = oa_i (:,:,:)
457      u_ice_b(:,:)     = u_ice(:,:)
458      v_ice_b(:,:)     = v_ice(:,:)
459      ! total concentration is needed for Lupkes parameterizations
460      at_i_b (:,:)     = at_i (:,:) 
461
462!!clem: output of initial state should be written here but it is impossible because
463!!      the ocean and ice are in the same file
464!!      CALL dia_wri_state( 'output.init' )
465      !
466   END SUBROUTINE ice_istate
467
468
469   SUBROUTINE ice_istate_init
470      !!-------------------------------------------------------------------
471      !!                   ***  ROUTINE ice_istate_init  ***
472      !!       
473      !! ** Purpose :   Definition of initial state of the ice
474      !!
475      !! ** Method  :   Read the namini namelist and check the parameter
476      !!              values called at the first timestep (nit000)
477      !!
478      !! ** input   :  Namelist namini
479      !!
480      !!-----------------------------------------------------------------------------
481      INTEGER ::   ios   ! Local integer output status for namelist read
482      INTEGER ::   ifpr, ierror
483      !
484      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
485      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
486      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
487      !
488      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
489         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
490         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
491         &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
492         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
493         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
494      !!-----------------------------------------------------------------------------
495      !
496      REWIND( numnam_ice_ref )              ! Namelist namini in reference namelist : Ice initial state
497      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
498901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namini in reference namelist' )
499      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
500      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
501902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namini in configuration namelist' )
502      IF(lwm) WRITE ( numoni, namini )
503      !
504      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
505      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
506      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
507      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd
508      !
509      IF(lwp) THEN                          ! control print
510         WRITE(numout,*)
511         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
512         WRITE(numout,*) '~~~~~~~~~~~~~~~'
513         WRITE(numout,*) '   Namelist namini:'
514         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini
515         WRITE(numout,*) '      ice initialization from a netcdf file            ln_iceini_file = ', ln_iceini_file
516         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
517         IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
518            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
519            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
520            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s
521            WRITE(numout,*) '      initial ice salinity  in the north-south         rn_smi_ini     = ', rn_smi_ini_n,rn_smi_ini_s
522            WRITE(numout,*) '      initial surf temperat in the north-south         rn_tsu_ini     = ', rn_tsu_ini_n,rn_tsu_ini_s
523            WRITE(numout,*) '      initial ice temperat  in the north-south         rn_tmi_ini     = ', rn_tmi_ini_n,rn_tmi_ini_s
524            WRITE(numout,*) '      initial snw temperat  in the north-south         rn_tms_ini     = ', rn_tms_ini_n,rn_tms_ini_s
525            WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
526            WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
527         ENDIF
528      ENDIF
529      !
530      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
531         !
532         ! set si structure
533         ALLOCATE( si(jpfldi), STAT=ierror )
534         IF( ierror > 0 ) THEN
535            CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' )   ;   RETURN
536         ENDIF
537         !
538         DO ifpr = 1, jpfldi
539            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
540            IF( slf_i(ifpr)%ln_tint )  ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
541         END DO
542         !
543         ! fill si with slf_i and control print
544         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
545         !
546      ENDIF
547      !
548      IF( .NOT.ln_pnd ) THEN
549         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
550         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
551         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
552      ENDIF
553      !
554   END SUBROUTINE ice_istate_init
555
556#else
557   !!----------------------------------------------------------------------
558   !!   Default option :         Empty module         NO SI3 sea-ice model
559   !!----------------------------------------------------------------------
560#endif
561
562   !!======================================================================
563END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.