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

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

changes in style - part5 - continue changing init routines

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