source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90 @ 8518

Last change on this file since 8518 was 8518, checked in by clem, 3 years ago

changes in style - part6 - commits of the day

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