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

source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceistate.F90 @ 13470

Last change on this file since 13470 was 13470, checked in by smasson, 4 years ago

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