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

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

source: NEMO/trunk/src/ICE/iceistate.F90 @ 9653

Last change on this file since 9653 was 9604, checked in by clem, 6 years ago

change history of the ice routines

File size: 29.2 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 icevar         ! ice_var_salprof
26   !
27   USE in_out_manager ! I/O manager
28   USE iom            ! I/O manager library
29   USE lib_mpp        ! MPP library
30   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
31   USE fldread        ! read input fields
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ice_istate        ! called by icestp.F90
37   PUBLIC   ice_istate_init   ! called by icestp.F90
38
39   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
40   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
41   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
42   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
43   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
44   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
45   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
46   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
47   !
48   !                             !! ** namelist (namini) **
49   LOGICAL  ::   ln_iceini        ! initialization or not
50   LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file
51   REAL(wp) ::   rn_thres_sst     ! threshold water temperature for initial sea ice
52   REAL(wp) ::   rn_hts_ini_n     ! initial snow thickness in the north
53   REAL(wp) ::   rn_hts_ini_s     ! initial snow thickness in the south
54   REAL(wp) ::   rn_hti_ini_n     ! initial ice thickness in the north
55   REAL(wp) ::   rn_hti_ini_s     ! initial ice thickness in the south
56   REAL(wp) ::   rn_ati_ini_n     ! initial leads area in the north
57   REAL(wp) ::   rn_ati_ini_s     ! initial leads area in the south
58   REAL(wp) ::   rn_smi_ini_n     ! initial salinity
59   REAL(wp) ::   rn_smi_ini_s     ! initial salinity
60   REAL(wp) ::   rn_tmi_ini_n     ! initial temperature
61   REAL(wp) ::   rn_tmi_ini_s     ! initial temperature
62   
63   !!----------------------------------------------------------------------
64   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
65   !! $Id: iceistate.F90 8378 2017-07-26 13:55:59Z clem $
66   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE ice_istate
71      !!-------------------------------------------------------------------
72      !!                    ***  ROUTINE ice_istate  ***
73      !!
74      !! ** Purpose :   defined the sea-ice initial state
75      !!
76      !! ** Method  :   This routine will put some ice where ocean
77      !!                is at the freezing point, then fill in ice
78      !!                state variables using prescribed initial
79      !!                values in the namelist           
80      !!
81      !! ** Steps   :   1) Set initial surface and basal temperatures
82      !!                2) Recompute or read sea ice state variables
83      !!                3) Fill in the ice thickness distribution using gaussian
84      !!                4) Fill in space-dependent arrays for state variables
85      !!                5) snow-ice mass computation
86      !!                6) store before fields
87      !!
88      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
89      !!              where there is no ice (clem: I do not know why, is it mandatory?)
90      !!--------------------------------------------------------------------
91      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
92      INTEGER  ::   i_hemis, i_fill, jl0   ! local integers
93      REAL(wp) ::   ztmelts, zdh
94      REAL(wp) ::   zarg, zV, zconv, zdv, zfac
95      INTEGER , DIMENSION(4)           ::   itest
96      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
97      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
98      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
99      REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
100      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini , za_i_ini                        !data by cattegories to fill
101      !--------------------------------------------------------------------
102
103      IF(lwp) WRITE(numout,*)
104      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
105      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
106
107      !--------------------------------------------------------------------
108      ! 1) Set surface and bottom temperatures to initial values
109      !--------------------------------------------------------------------
110      !
111      ! init surface temperature
112      DO jl = 1, jpl
113         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface
114         cnd_ice(:,:,jl) = 0._wp               ! initialisation of the effective conductivity at the top of ice/snow (Jules coupling)
115      END DO
116      !
117      ! init basal temperature (considered at freezing point)   [Kelvin]
118      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
119      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
120
121      IF( ln_iceini ) THEN
122         !-----------------------------------------------------------
123         ! 2) Compute or read sea ice variables ===> single category
124         !-----------------------------------------------------------
125         !
126         !                             !---------------!
127         IF( ln_iceini_file )THEN      ! Read a file   !
128            !                          !---------------!
129            !
130            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
131            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
132            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
133            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
134            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
135            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
136            !
137            WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 
138            ELSEWHERE                       ; zswitch(:,:) = 0._wp
139            END WHERE
140            zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)
141            !
142            !                          !---------------!
143         ELSE                          ! Read namelist !
144            !                          !---------------!
145            ! no ice if sst <= t-freez + ttest
146            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
147            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
148            END WHERE
149            !
150            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
151            WHERE( ff_t(:,:) >= 0._wp )
152               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
153               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
154               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
155               zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
156               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
157               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
158            ELSEWHERE
159               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
160               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
161               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
162               zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
163               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
164               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
165            END WHERE
166            zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)
167            !
168         ENDIF
169         
170         !------------------------------------------------------------------
171         ! 3) Distribute ice concentration and thickness into the categories
172         !------------------------------------------------------------------
173         ! a gaussian distribution for ice concentration is used
174         ! then we check whether the distribution fullfills
175         ! volume and area conservation, positivity and ice categories bounds
176         zh_i_ini(:,:,:) = 0._wp 
177         za_i_ini(:,:,:) = 0._wp
178         !
179         DO jj = 1, jpj
180            DO ji = 1, jpi
181               !
182               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
183
184                  ! find which category (jl0) the input ice thickness falls into
185                  jl0 = jpl
186                  DO jl = 1, jpl
187                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
188                        jl0 = jl
189                        CYCLE
190                     ENDIF
191                  END DO
192                  !
193                  itest(:) = 0
194                  i_fill   = jpl + 1                                            !------------------------------------
195                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
196                     !                                                          !------------------------------------
197                     i_fill = i_fill - 1
198                     !
199                     zh_i_ini(ji,jj,:) = 0._wp 
200                     za_i_ini(ji,jj,:) = 0._wp
201                     itest(:) = 0
202                     !
203                     IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
204                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
205                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
206                     ELSE                         !-- case ice is thicker: fill categories >1
207                        ! thickness
208                        DO jl = 1, i_fill-1
209                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
210                        END DO
211                        !
212                        ! concentration
213                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
214                        DO jl = 1, i_fill - 1
215                           IF( jl /= jl0 )THEN
216                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
217                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
218                           ENDIF
219                        END DO
220
221                        ! last category
222                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
223                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
224                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
225
226                        ! correction if concentration of upper cat is greater than lower cat
227                        !   (it should be a gaussian around jl0 but sometimes it is not)
228                        IF ( jl0 /= jpl ) THEN
229                           DO jl = jpl, jl0+1, -1
230                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
231                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
232                                 zh_i_ini(ji,jj,jl    ) = 0._wp
233                                 za_i_ini(ji,jj,jl    ) = 0._wp
234                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
235                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
236                              END IF
237                           ENDDO
238                        ENDIF
239                        !
240                     ENDIF
241                     !
242                     ! Compatibility tests
243                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation
244                     IF ( zconv < epsi06 ) itest(1) = 1
245                     !
246                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation
247                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
248                     IF ( zconv < epsi06 ) itest(2) = 1
249                     !
250                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ?
251                     !
252                     itest(4) = 1
253                     DO jl = 1, i_fill
254                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations
255                     END DO
256                     !                                                          !----------------------------
257                  END DO                                                        ! end iteration on categories
258                  !                                                             !----------------------------
259                  IF( lwp .AND. SUM(itest) /= 4 ) THEN
260                     WRITE(numout,*)
261                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
262                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
263                     WRITE(numout,*)
264                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
265                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
266                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
267                  ENDIF
268                  !
269               ENDIF
270               !
271            END DO   
272         END DO   
273
274         !---------------------------------------------------------------------
275         ! 4) Fill in sea ice arrays
276         !---------------------------------------------------------------------
277         !
278         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
279         DO jl = 1, jpl ! loop over categories
280            DO jj = 1, jpj
281               DO ji = 1, jpi
282                  a_i(ji,jj,jl)  = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
283                  h_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
284                  s_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
285                  o_i(ji,jj,jl)  = 0._wp                                                     ! age (0 day)
286                  t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
287                  !
288                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
289                    h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
290                  ELSE
291                    h_s(ji,jj,jl)= 0._wp
292                  ENDIF
293                  !
294                  ! This case below should not be used if (h_s/h_i) is ok in namelist
295                  ! In case snow load is in excess that would lead to transformation from snow to ice
296                  ! Then, transfer the snow excess into the ice (different from icethd_dh)
297                  zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 
298                  ! recompute h_i, h_s avoiding out of bounds values
299                  h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
300                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
301                  !
302                  ! ice volume, salt content, age content
303                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
304                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
305                  sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
306                  oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
307               END DO
308            END DO
309         END DO
310         !
311         IF( nn_icesal /= 2 )  THEN         ! for constant salinity in time
312            CALL ice_var_salprof
313            sv_i = s_i * v_i
314         ENDIF
315         
316         ! Snow temperature and heat content
317         DO jk = 1, nlay_s
318            DO jl = 1, jpl ! loop over categories
319               DO jj = 1, jpj
320                  DO ji = 1, jpi
321                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
322                     ! Snow energy of melting
323                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
324                     !
325                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
326                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
327                  END DO
328               END DO
329            END DO
330         END DO
331         !
332         ! Ice salinity, temperature and heat content
333         DO jk = 1, nlay_i
334            DO jl = 1, jpl ! loop over categories
335               DO jj = 1, jpj
336                  DO ji = 1, jpi
337                     t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
338                     sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
339                     ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
340                     !
341                     ! heat content per unit volume
342                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) )         &
343                        &             + lfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   &
344                        &             - rcp  * ( ztmelts - rt0 ) )
345                     !
346                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
347                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
348                  END DO
349               END DO
350            END DO
351         END DO
352         !
353         tn_ice (:,:,:) = t_su (:,:,:)
354         t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu
355
356         ! Melt pond volume and fraction
357         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp
358         ELSE                                     ;   zfac = 0._wp
359         ENDIF
360         DO jl = 1, jpl
361            a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac
362            h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac
363         END DO
364         a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 
365         v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:)
366         !
367      ELSE ! if ln_iceini=false
368         a_i  (:,:,:) = 0._wp
369         v_i  (:,:,:) = 0._wp
370         v_s  (:,:,:) = 0._wp
371         sv_i (:,:,:) = 0._wp
372         oa_i (:,:,:) = 0._wp
373         h_i  (:,:,:) = 0._wp
374         h_s  (:,:,:) = 0._wp
375         s_i  (:,:,:) = 0._wp
376         o_i  (:,:,:) = 0._wp
377         !
378         e_i(:,:,:,:) = 0._wp
379         e_s(:,:,:,:) = 0._wp
380         !
381         DO jl = 1, jpl
382            DO jk = 1, nlay_i
383               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
384            END DO
385            DO jk = 1, nlay_s
386               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
387            END DO
388         END DO
389
390         tn_ice (:,:,:) = t_i (:,:,1,:)
391         t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu
392         
393         a_ip(:,:,:)      = 0._wp
394         v_ip(:,:,:)      = 0._wp
395         a_ip_frac(:,:,:) = 0._wp
396         h_ip     (:,:,:) = 0._wp
397         !
398      ENDIF ! ln_iceini
399      !
400      at_i (:,:) = 0.0_wp
401      DO jl = 1, jpl
402         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
403      END DO
404      !
405      ! --- set ice velocities --- !
406      u_ice (:,:) = 0._wp
407      v_ice (:,:) = 0._wp
408      !
409      !----------------------------------------------
410      ! 5) Snow-ice mass (case ice is fully embedded)
411      !----------------------------------------------
412      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhosn * v_s(:,:,:) + rhoic * v_i(:,:,:), dim=3  )   ! snow+ice mass
413      snwice_mass_b(:,:) = snwice_mass(:,:)
414      !
415      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
416         !
417         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
418         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
419         !
420         IF( .NOT.ln_linssh ) THEN
421            !
422            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
423            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
424            !
425            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
426               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
427               e3t_b(:,:,jk) = e3t_n(:,:,jk)
428               e3t_a(:,:,jk) = e3t_n(:,:,jk)
429            END DO
430            !
431            ! Reconstruction of all vertical scale factors at now and before time-steps
432            ! =========================================================================
433            ! Horizontal scale factor interpolations
434            ! --------------------------------------
435            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
436            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
437            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
438            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
439            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
440            ! Vertical scale factor interpolations
441            ! ------------------------------------
442            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
443            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
444            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
445            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
446            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
447            ! t- and w- points depth
448            ! ----------------------
449            !!gm not sure of that....
450            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
451            gdepw_n(:,:,1) = 0.0_wp
452            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
453            DO jk = 2, jpk
454               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
455               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
456               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
457            END DO
458         ENDIF
459      ENDIF
460     
461      !------------------------------------
462      ! 6) store fields at before time-step
463      !------------------------------------
464      ! it is only necessary for the 1st interpolation by Agrif
465      a_i_b  (:,:,:)   = a_i  (:,:,:)
466      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
467      v_i_b  (:,:,:)   = v_i  (:,:,:)
468      v_s_b  (:,:,:)   = v_s  (:,:,:)
469      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
470      sv_i_b (:,:,:)   = sv_i (:,:,:)
471      oa_i_b (:,:,:)   = oa_i (:,:,:)
472      u_ice_b(:,:)     = u_ice(:,:)
473      v_ice_b(:,:)     = v_ice(:,:)
474
475!!clem: output of initial state should be written here but it is impossible because
476!!      the ocean and ice are in the same file
477!!      CALL dia_wri_state( 'output.init', nit000 )
478      !
479   END SUBROUTINE ice_istate
480
481
482   SUBROUTINE ice_istate_init
483      !!-------------------------------------------------------------------
484      !!                   ***  ROUTINE ice_istate_init  ***
485      !!       
486      !! ** Purpose :   Definition of initial state of the ice
487      !!
488      !! ** Method  :   Read the namini namelist and check the parameter
489      !!              values called at the first timestep (nit000)
490      !!
491      !! ** input   :  Namelist namini
492      !!
493      !!-----------------------------------------------------------------------------
494      INTEGER ::   ji, jj
495      INTEGER ::   ios, ierr, inum_ice   ! Local integer output status for namelist read
496      INTEGER ::   ifpr, ierror
497      !
498      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
499      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
500      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
501      !
502      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
503         &             rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
504         &             rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
505         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
506      !!-----------------------------------------------------------------------------
507      !
508      REWIND( numnam_ice_ref )              ! Namelist namini in reference namelist : Ice initial state
509      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
510901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namini in reference namelist', lwp )
511      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
512      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
513902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namini in configuration namelist', lwp )
514      IF(lwm) WRITE ( numoni, namini )
515      !
516      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
517      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
518      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
519      !
520      IF(lwp) THEN                          ! control print
521         WRITE(numout,*)
522         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
523         WRITE(numout,*) '~~~~~~~~~~~~~~~'
524         WRITE(numout,*) '   Namelist namini:'
525         WRITE(numout,*) '      initialization with ice (T) or not (F)                 ln_iceini       = ', ln_iceini
526         WRITE(numout,*) '      ice initialization from a netcdf file                  ln_iceini_file  = ', ln_iceini_file
527         WRITE(numout,*) '      max delta ocean temp. above Tfreeze with initial ice   rn_thres_sst    = ', rn_thres_sst
528         WRITE(numout,*) '      initial snow thickness in the north                    rn_hts_ini_n    = ', rn_hts_ini_n
529         WRITE(numout,*) '      initial snow thickness in the south                    rn_hts_ini_s    = ', rn_hts_ini_s 
530         WRITE(numout,*) '      initial ice thickness  in the north                    rn_hti_ini_n    = ', rn_hti_ini_n
531         WRITE(numout,*) '      initial ice thickness  in the south                    rn_hti_ini_s    = ', rn_hti_ini_s
532         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_n    = ', rn_ati_ini_n
533         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_s    = ', rn_ati_ini_s
534         WRITE(numout,*) '      initial  ice salinity  in the north                    rn_smi_ini_n    = ', rn_smi_ini_n
535         WRITE(numout,*) '      initial  ice salinity  in the south                    rn_smi_ini_s    = ', rn_smi_ini_s
536         WRITE(numout,*) '      initial  ice/snw temp  in the north                    rn_tmi_ini_n    = ', rn_tmi_ini_n
537         WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s    = ', rn_tmi_ini_s
538      ENDIF
539      !
540      IF( ln_iceini_file ) 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_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            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', 'ice istate ini', 'numnam_ice' )
555         !
556         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
557         !
558      ENDIF
559      !
560   END SUBROUTINE ice_istate_init
561
562#else
563   !!----------------------------------------------------------------------
564   !!   Default option :         Empty module         NO SI3 sea-ice model
565   !!----------------------------------------------------------------------
566#endif
567
568   !!======================================================================
569END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.