New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limistate.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 7881

Last change on this file since 7881 was 7761, checked in by clem, 7 years ago

make AGRIF and LIM3 fully compatible

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