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

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 26.6 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   !! * Substitutions
64#  include "do_loop_substitute.h90"
65   !!----------------------------------------------------------------------
66   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
67   !! $Id$
68   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
69   !!----------------------------------------------------------------------
70CONTAINS
71
72   SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa )
73      !!-------------------------------------------------------------------
74      !!                    ***  ROUTINE ice_istate  ***
75      !!
76      !! ** Purpose :   defined the sea-ice initial state
77      !!
78      !! ** Method  :   This routine will put some ice where ocean
79      !!                is at the freezing point, then fill in ice
80      !!                state variables using prescribed initial
81      !!                values in the namelist           
82      !!
83      !! ** Steps   :   1) Set initial surface and basal temperatures
84      !!                2) Recompute or read sea ice state variables
85      !!                3) Fill in the ice thickness distribution using gaussian
86      !!                4) Fill in space-dependent arrays for state variables
87      !!                5) snow-ice mass computation
88      !!                6) store before fields
89      !!
90      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
91      !!              where there is no ice
92      !!--------------------------------------------------------------------
93      INTEGER, INTENT(in) :: kt            ! time step
94      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
95      !
96      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
97      REAL(wp) ::   ztmelts
98      INTEGER , DIMENSION(4)           ::   itest
99      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
100      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
101      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
102      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
103      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file
104      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
105      !!
106      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
107      !--------------------------------------------------------------------
108
109      IF(lwp) WRITE(numout,*)
110      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
111      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
112
113      !---------------------------
114      ! 1) 1st init. of the fields
115      !---------------------------
116      !
117      ! basal temperature (considered at freezing point)   [Kelvin]
118      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
119      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
120      !
121      ! surface temperature and conductivity
122      DO jl = 1, jpl
123         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface
124         cnd_ice(:,:,jl) = 0._wp               ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
125      END DO
126      !
127      ! ice and snw temperatures
128      DO jl = 1, jpl
129         DO jk = 1, nlay_i
130            t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
131         END DO
132         DO jk = 1, nlay_s
133            t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
134         END DO
135      END DO
136      !
137      ! specific temperatures for coupled runs
138      tn_ice (:,:,:) = t_i (:,:,1,:)
139      t1_ice (:,:,:) = t_i (:,:,1,:)
140
141      ! heat contents
142      e_i (:,:,:,:) = 0._wp
143      e_s (:,:,:,:) = 0._wp
144     
145      ! general fields
146      a_i (:,:,:) = 0._wp
147      v_i (:,:,:) = 0._wp
148      v_s (:,:,:) = 0._wp
149      sv_i(:,:,:) = 0._wp
150      oa_i(:,:,:) = 0._wp
151      !
152      h_i (:,:,:) = 0._wp
153      h_s (:,:,:) = 0._wp
154      s_i (:,:,:) = 0._wp
155      o_i (:,:,:) = 0._wp
156      !
157      ! melt ponds
158      a_ip     (:,:,:) = 0._wp
159      v_ip     (:,:,:) = 0._wp
160      a_ip_frac(:,:,:) = 0._wp
161      h_ip     (:,:,:) = 0._wp
162      !
163      ! ice velocities
164      u_ice (:,:) = 0._wp
165      v_ice (:,:) = 0._wp
166      !
167      !------------------------------------------------------------------------
168      ! 2) overwrite some of the fields with namelist parameters or netcdf file
169      !------------------------------------------------------------------------
170      IF( ln_iceini ) THEN
171         !                             !---------------!
172         IF( ln_iceini_file )THEN      ! Read a file   !
173            !                          !---------------!
174            WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
175            ELSEWHERE                     ;   zswitch(:,:) = 0._wp
176            END WHERE
177            !
178            CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
179            !
180            ! -- mandatory fields -- !
181            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
182            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
183            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
184
185            ! -- optional fields -- !
186            !    if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
187            !
188            ! ice salinity
189            IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
190               &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
191            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
192            !
193            ! ice temperature
194            IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) &
195               &     si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
196            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
197            !
198            ! surface temperature => set to ice temperature if it exists
199            IF    ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
200                     si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
201            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
202                     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
203            ENDIF
204            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
205            !
206            ! snow temperature => set to ice temperature if it exists
207            IF    ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
208                     si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
209            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
210                     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
211            ENDIF
212            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
213            !
214            ! pond concentration
215            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
216               &     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.
217               &                              * si(jp_ati)%fnow(:,:,1) 
218            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
219            !
220            ! pond depth
221            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
222               &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
223            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,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_2D_11_11
274            IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
275               npti         = npti  + 1
276               nptidx(npti) = (jj - 1) * jpi + ji
277            ENDIF
278         END_2D
279
280         ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
281         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
282         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
283         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
284         CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
285         CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
286         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
287         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
288         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
289         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
290
291         ! allocate temporary arrays
292         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
293            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
294         
295         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
296         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
297            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   &
298            &              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), &
299            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d )
300
301         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
302         DO jl = 1, jpl
303            zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
304            zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
305         END DO
306         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
307         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
308         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
309         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
310         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
311         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
312         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
313         CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
314         CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
315
316         ! deallocate temporary arrays
317         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
318            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
319
320         ! calculate extensive and intensive variables
321         CALL ice_var_salprof ! for sz_i
322         DO jl = 1, jpl
323            DO_2D_11_11
324               v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
325               v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
326               sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
327            END_2D
328         END DO
329         !
330         DO jl = 1, jpl
331            DO_3D_11_11( 1, nlay_s )
332               t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
333               e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
334                  &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
335            END_3D
336         END DO
337         !
338         DO jl = 1, jpl
339            DO_3D_11_11( 1, nlay_i )
340               t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
341               ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
342               e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
343                  &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
344                  &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
345                  &                       - rcp   * ( ztmelts - rt0 ) )
346            END_3D
347         END DO
348
349         ! Melt ponds
350         WHERE( a_i > epsi10 )
351            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
352         ELSEWHERE
353            a_ip_frac(:,:,:) = 0._wp
354         END WHERE
355         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
356         
357         ! specific temperatures for coupled runs
358         tn_ice(:,:,:) = t_su(:,:,:)
359         t1_ice(:,:,:) = t_i (:,:,1,:)
360         !
361      ENDIF ! ln_iceini
362      !
363      at_i(:,:) = SUM( a_i, dim=3 )
364      !
365      !----------------------------------------------
366      ! 3) Snow-ice mass (case ice is fully embedded)
367      !----------------------------------------------
368      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass
369      snwice_mass_b(:,:) = snwice_mass(:,:)
370      !
371      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
372         !
373         ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rau0
374         ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rau0
375         !
376         IF( .NOT.ln_linssh ) THEN
377            !
378            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:)
379            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
380            !
381            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
382               e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:)
383               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
384               e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm)
385            END DO
386            !
387            ! Reconstruction of all vertical scale factors at now and before time-steps
388            ! =========================================================================
389            ! Horizontal scale factor interpolations
390            ! --------------------------------------
391            CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
392            CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
393            CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
394            CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
395            CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
396            ! Vertical scale factor interpolations
397            ! ------------------------------------
398            CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
399            CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
400            CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
401            CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
402            CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
403            ! t- and w- points depth
404            ! ----------------------
405            !!gm not sure of that....
406            gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
407            gdepw(:,:,1,Kmm) = 0.0_wp
408            gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
409            DO jk = 2, jpk
410               gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk  ,Kmm)
411               gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
412               gde3w(:,:,jk) = gdept(:,:,jk  ,Kmm) - ssh (:,:,Kmm)
413            END DO
414         ENDIF
415      ENDIF
416     
417      !------------------------------------
418      ! 4) store fields at before time-step
419      !------------------------------------
420      ! it is only necessary for the 1st interpolation by Agrif
421      a_i_b  (:,:,:)   = a_i  (:,:,:)
422      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
423      v_i_b  (:,:,:)   = v_i  (:,:,:)
424      v_s_b  (:,:,:)   = v_s  (:,:,:)
425      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
426      sv_i_b (:,:,:)   = sv_i (:,:,:)
427      oa_i_b (:,:,:)   = oa_i (:,:,:)
428      u_ice_b(:,:)     = u_ice(:,:)
429      v_ice_b(:,:)     = v_ice(:,:)
430      ! total concentration is needed for Lupkes parameterizations
431      at_i_b (:,:)     = at_i (:,:) 
432
433!!clem: output of initial state should be written here but it is impossible because
434!!      the ocean and ice are in the same file
435!!      CALL dia_wri_state( 'output.init' )
436      !
437   END SUBROUTINE ice_istate
438
439
440   SUBROUTINE ice_istate_init
441      !!-------------------------------------------------------------------
442      !!                   ***  ROUTINE ice_istate_init  ***
443      !!       
444      !! ** Purpose :   Definition of initial state of the ice
445      !!
446      !! ** Method  :   Read the namini namelist and check the parameter
447      !!              values called at the first timestep (nit000)
448      !!
449      !! ** input   :  Namelist namini
450      !!
451      !!-----------------------------------------------------------------------------
452      INTEGER ::   ios   ! Local integer output status for namelist read
453      INTEGER ::   ifpr, ierror
454      !
455      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
456      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
457      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
458      !
459      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
460         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
461         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
462         &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
463         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
464         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
465      !!-----------------------------------------------------------------------------
466      !
467      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
468901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namini in reference namelist' )
469      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
470902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namini in configuration namelist' )
471      IF(lwm) WRITE ( numoni, namini )
472      !
473      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
474      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
475      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
476      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd
477      !
478      IF(lwp) THEN                          ! control print
479         WRITE(numout,*)
480         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
481         WRITE(numout,*) '~~~~~~~~~~~~~~~'
482         WRITE(numout,*) '   Namelist namini:'
483         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini
484         WRITE(numout,*) '      ice initialization from a netcdf file            ln_iceini_file = ', ln_iceini_file
485         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
486         IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
487            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
488            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
489            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s
490            WRITE(numout,*) '      initial ice salinity  in the north-south         rn_smi_ini     = ', rn_smi_ini_n,rn_smi_ini_s
491            WRITE(numout,*) '      initial surf temperat in the north-south         rn_tsu_ini     = ', rn_tsu_ini_n,rn_tsu_ini_s
492            WRITE(numout,*) '      initial ice temperat  in the north-south         rn_tmi_ini     = ', rn_tmi_ini_n,rn_tmi_ini_s
493            WRITE(numout,*) '      initial snw temperat  in the north-south         rn_tms_ini     = ', rn_tms_ini_n,rn_tms_ini_s
494            WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
495            WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
496         ENDIF
497      ENDIF
498      !
499      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
500         !
501         ! set si structure
502         ALLOCATE( si(jpfldi), STAT=ierror )
503         IF( ierror > 0 ) THEN
504            CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' )   ;   RETURN
505         ENDIF
506         !
507         DO ifpr = 1, jpfldi
508            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
509            IF( slf_i(ifpr)%ln_tint )  ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
510         END DO
511         !
512         ! fill si with slf_i and control print
513         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
514         !
515      ENDIF
516      !
517      IF( .NOT.ln_pnd ) THEN
518         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
519         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
520         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
521      ENDIF
522      !
523   END SUBROUTINE ice_istate_init
524
525#else
526   !!----------------------------------------------------------------------
527   !!   Default option :         Empty module         NO SI3 sea-ice model
528   !!----------------------------------------------------------------------
529#endif
530
531   !!======================================================================
532END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.