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

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

changes in style - part5 - very nearly finished

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