source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 8291

Last change on this file since 8291 was 8291, checked in by clem, 4 years ago

correct last commit changes. This is the start version for the modifications of the ice model. All sette tests passed except restartability for agrif and sas (though sas_biper is restartable). Configurations CREG, SPITZ, FER & SASBIPER are running.

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