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

Last change on this file since 13180 was 13180, checked in by smueller, 4 years ago

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