New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
iceistate.F90 in NEMO/branches/2020/ticket2487/src/ICE – NEMO

source: NEMO/branches/2020/ticket2487/src/ICE/iceistate.F90 @ 15410

Last change on this file since 15410 was 15410, checked in by smueller, 12 months ago

Synchronizing with /NEMO/releases/r4.0/r4.0-HEAD@15405 (ticket #2487)

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