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/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceistate.F90 @ 15451

Last change on this file since 15451 was 15451, checked in by clem, 8 months ago

4.0-HEAD only: change vertical scale factors for embedded sea ice as proposed in ticket #2487.

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