source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/ICE/iceistate.F90 @ 11021

Last change on this file since 11021 was 10998, checked in by davestorkey, 22 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : ICE changes. Passess non-AGRIF sette tests apart from ORCA2_SAS_ICE which doesn't seem to work for me at the moment.

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