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

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90 @ 12734

Last change on this file since 12734 was 12728, checked in by clem, 4 years ago

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