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/UKMO/NEMO_4.0_mirror/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/iceistate.F90 @ 10888

Last change on this file since 10888 was 10888, checked in by davestorkey, 5 years ago

branches/UKMO/NEMO_4.0_mirror : clear SVN keywords

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