source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90 @ 8498

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

changes in style - part2 -

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