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

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

last routine names to be changed

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