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 branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 9 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File size: 23.8 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
[4161]7   !!            4.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
[2528]25   USE dom_ice          ! sea-ice domain
26   USE in_out_manager   ! I/O manager
[2715]27   USE lib_mpp          ! MPP library
[4161]28   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[3294]29   USE wrk_nemo         ! work arrays
[825]30
31   IMPLICIT NONE
32   PRIVATE
33
[2528]34   PUBLIC   lim_istate      ! routine called by lim_init.F90
[825]35
[4166]36   !                          !!** init namelist (namiceini) **
[5443]37   REAL(wp) ::   rn_thres_sst   ! threshold water temperature for initial sea ice
38   REAL(wp) ::   rn_hts_ini_n   ! initial snow thickness in the north
39   REAL(wp) ::   rn_hts_ini_s   ! initial snow thickness in the south
40   REAL(wp) ::   rn_hti_ini_n   ! initial ice thickness in the north
41   REAL(wp) ::   rn_hti_ini_s   ! initial ice thickness in the south
42   REAL(wp) ::   rn_ati_ini_n   ! initial leads area in the north
43   REAL(wp) ::   rn_ati_ini_s   ! initial leads area in the south
44   REAL(wp) ::   rn_smi_ini_n   ! initial salinity
45   REAL(wp) ::   rn_smi_ini_s   ! initial salinity
46   REAL(wp) ::   rn_tmi_ini_n   ! initial temperature
47   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature
[4166]48
[5443]49   LOGICAL  ::  ln_iceini    ! initialization or not
[825]50   !!----------------------------------------------------------------------
[4161]51   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
[5234]52   !! $Id$
[4161]53   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[825]54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE lim_istate
58      !!-------------------------------------------------------------------
59      !!                    ***  ROUTINE lim_istate  ***
60      !!
61      !! ** Purpose :   defined the sea-ice initial state
62      !!
[4161]63      !! ** Method  :   
64      !!                This routine will put some ice where ocean
65      !!                is at the freezing point, then fill in ice
66      !!                state variables using prescribed initial
67      !!                values in the namelist           
68      !!
69      !! ** Steps   :   
70      !!                1) Read namelist
71      !!                2) Basal temperature; ice and hemisphere masks
72      !!                3) Fill in the ice thickness distribution using gaussian
73      !!                4) Fill in space-dependent arrays for state variables
74      !!                5) Diagnostic arrays
75      !!                6) Lateral boundary conditions
76      !!
[4333]77      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
[4990]78      !!              where there is no ice (clem: I do not know why, is it mandatory?)
[4333]79      !!
[4161]80      !! History :
81      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
82      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
83      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
84      !!--------------------------------------------------------------------
85
86      !! * Local variables
87      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
[5443]88      REAL(wp)   :: ztmelts, zdh
[4161]89      INTEGER    :: i_hemis, i_fill, jl0 
90      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
[4688]91      REAL(wp), POINTER, DIMENSION(:)     :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini
92      REAL(wp), POINTER, DIMENSION(:,:)   :: zh_i_ini, za_i_ini, zv_i_ini
93      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator
[4161]94      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index
[825]95      !--------------------------------------------------------------------
[921]96
[4688]97      CALL wrk_alloc( jpi, jpj, zswitch )
[4161]98      CALL wrk_alloc( jpi, jpj, zhemis )
[4688]99      CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
100      CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
[2715]101
[4161]102      IF(lwp) WRITE(numout,*)
103      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
104      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
105
[825]106      !--------------------------------------------------------------------
[4161]107      ! 1) Read namelist
[825]108      !--------------------------------------------------------------------
109
110      CALL lim_istate_init     !  reading the initials parameters of the ice
111
[4688]112      ! surface temperature
113      DO jl = 1, jpl ! loop over categories
[5443]114         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
115         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
[4688]116      END DO
117
[4990]118      ! basal temperature (considered at freezing point)
[5443]119      t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1) 
[4990]120
[5443]121      IF( ln_iceini ) THEN
[4688]122
[825]123      !--------------------------------------------------------------------
[4161]124      ! 2) Basal temperature, ice mask and hemispheric index
[825]125      !--------------------------------------------------------------------
[4990]126
127      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
[1037]128         DO ji = 1, jpi
[5443]129            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN
130               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
[4765]131            ELSE                                                                                   
[5443]132               zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
[1037]133            ENDIF
134         END DO
135      END DO
136
137
[4161]138      ! Hemispheric index
[825]139      DO jj = 1, jpj
140         DO ji = 1, jpi
[2528]141            IF( fcor(ji,jj) >= 0._wp ) THEN   
[4161]142               zhemis(ji,jj) = 1 ! Northern hemisphere
143            ELSE
144               zhemis(ji,jj) = 2 ! Southern hemisphere
145            ENDIF
146         END DO
147      END DO
[825]148
[4161]149      !--------------------------------------------------------------------
150      ! 3) Initialization of sea ice state variables
151      !--------------------------------------------------------------------
[825]152
[4161]153      !-----------------------------
154      ! 3.1) Hemisphere-dependent arrays
155      !-----------------------------
[4688]156      ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
[5443]157      zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness
158      zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth
159      zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration
160      zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity
161      ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow)
[825]162
[4688]163      zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume
164
[4161]165      !---------------------------------------------------------------------
166      ! 3.2) Distribute ice concentration and thickness into the categories
167      !---------------------------------------------------------------------
168      ! a gaussian distribution for ice concentration is used
169      ! then we check whether the distribution fullfills
170      ! volume and area conservation, positivity and ice categories bounds
171      DO i_hemis = 1, 2
[825]172
[4161]173      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
[825]174
[4161]175      ! note for the great nemo engineers:
176      ! only very few of the WRITE statements are necessary for the reference version
177      ! they were one day useful, but now i personally doubt of their
178      ! potential for bringing anything useful
[825]179
[4161]180      DO i_fill = jpl, 1, -1
181         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
182            !----------------------------
183            ! fill the i_fill categories
184            !----------------------------
185            ! *** 1 category to fill
186            IF ( i_fill .EQ. 1 ) THEN
[4688]187               zh_i_ini(1,i_hemis)       = zht_i_ini(i_hemis)
188               za_i_ini(1,i_hemis)       = zat_i_ini(i_hemis)
189               zh_i_ini(2:jpl,i_hemis)   = 0._wp
190               za_i_ini(2:jpl,i_hemis)   = 0._wp
[4161]191            ELSE
[825]192
[4688]193               ! *** >1 categores to fill
194               !--- Ice thicknesses in the i_fill - 1 first categories
[4161]195               DO jl = 1, i_fill - 1
[5443]196                  zh_i_ini(jl,i_hemis) = hi_mean(jl)
[4161]197               END DO
[4688]198               
199               !--- jl0: most likely index where cc will be maximum
[825]200               DO jl = 1, jpl
[5443]201                  IF ( ( zht_i_ini(i_hemis) >  hi_max(jl-1) ) .AND. &
202                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN
[4161]203                     jl0 = jl
204                  ENDIF
205               END DO
206               jl0 = MIN(jl0, i_fill)
[4688]207               
208               !--- Concentrations
[4161]209               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl))
210               DO jl = 1, i_fill - 1
211                  IF ( jl .NE. jl0 ) THEN
[4688]212                     zsigma               = 0.5 * zht_i_ini(i_hemis)
213                     zarg                 = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma
[4161]214                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2)
215                  ENDIF
[4688]216               END DO
217               
[4161]218               zA = 0. ! sum of the areas in the jpl categories
219               DO jl = 1, i_fill - 1
220                 zA = zA + za_i_ini(jl,i_hemis)
221               END DO
222               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category
223               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
224         
[4688]225               !--- Ice thickness in the last category
[4161]226               zV = 0. ! sum of the volumes of the N-1 categories
227               DO jl = 1, i_fill - 1
[4688]228                  zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis)
[4161]229               END DO
[4688]230               zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 
231               IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
[825]232
[4688]233               !--- volumes
234               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis)
[4161]235               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
[921]236
[4161]237            ENDIF ! i_fill
[825]238
[4161]239            !---------------------
240            ! Compatibility tests
241            !---------------------
242            ! Test 1: area conservation
243            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons )
244            IF ( zconv .LT. 1.0e-6 ) THEN
245               ztest_1 = 1
246            ELSE 
247              ! this write is useful
248              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 
249               ztest_1 = 0
250            ENDIF
[825]251
[4161]252            ! Test 2: volume conservation
253            zV_cons = SUM(zv_i_ini(:,i_hemis))
254            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons)
[825]255
[4161]256            IF ( zconv .LT. 1.0e-6 ) THEN
257               ztest_2 = 1
258            ELSE
259              ! this write is useful
260              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
261                            ' zvt_i_ini = ', zvt_i_ini(i_hemis)
262               ztest_2 = 0
263            ENDIF
[825]264
[4161]265            ! Test 3: thickness of the last category is in-bounds ?
[5443]266            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN
[4161]267               ztest_3 = 1
268            ELSE
269               ! this write is useful
[4688]270               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', &
271               zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
[4161]272               ztest_3 = 0
273            ENDIF
[825]274
[4161]275            ! Test 4: positivity of ice concentrations
276            ztest_4 = 1
277            DO jl = 1, jpl
278               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
279                  ! this write is useful
280                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
281                  ztest_4 = 0
[825]282               ENDIF
[4161]283            END DO
[825]284
[4161]285         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
286 
287         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
[825]288
[4161]289      END DO ! i_fill
[825]290
[4161]291      IF(lwp) THEN
[4908]292         WRITE(numout,*) ' ztests : ', ztests
[4161]293         IF ( ztests .NE. 4 ) THEN
294            WRITE(numout,*)
[4908]295            WRITE(numout,*) ' !!!! ALERT                  !!! '
296            WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
[4161]297            WRITE(numout,*)
[4908]298            WRITE(numout,*) ' *** ztests is not equal to 4 '
299            WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
300            WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis)
301            WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis)
[4161]302         ENDIF ! ztests .NE. 4
303      ENDIF
304     
305      END DO ! i_hemis
[825]306
[4161]307      !---------------------------------------------------------------------
308      ! 3.3) Space-dependent arrays for ice state variables
309      !---------------------------------------------------------------------
[825]310
[4333]311      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
[4161]312      DO jl = 1, jpl ! loop over categories
313         DO jj = 1, jpj
314            DO ji = 1, jpi
[4688]315               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
[5443]316               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness
[4688]317               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth
[5443]318               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
319               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day)
320               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
[4333]321
322               ! This case below should not be used if (ht_s/ht_i) is ok in namelist
323               ! In case snow load is in excess that would lead to transformation from snow to ice
324               ! Then, transfer the snow excess into the ice (different from limthd_dh)
325               zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
326               ! recompute ht_i, ht_s avoiding out of bounds values
327               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
[5443]328               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
[4333]329
330               ! ice volume, salt content, age content
[4161]331               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
332               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
333               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
334               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
[5443]335            END DO
336         END DO
337      END DO
[825]338
[4161]339      ! Snow temperature and heat content
340      DO jk = 1, nlay_s
341         DO jl = 1, jpl ! loop over categories
342            DO jj = 1, jpj
343               DO ji = 1, jpi
[5443]344                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
[4161]345                   ! Snow energy of melting
[5443]346                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
[825]347
[5443]348                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
349                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
350               END DO
351            END DO
352         END DO
353      END DO
354
[4161]355      ! Ice salinity, temperature and heat content
356      DO jk = 1, nlay_i
357         DO jl = 1, jpl ! loop over categories
358            DO jj = 1, jpj
359               DO ji = 1, jpi
[5443]360                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
361                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
362                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
[825]363
[4161]364                   ! heat content per unit volume
[4688]365                   e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
[5443]366                      +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
367                      -   rcp     * ( ztmelts - rt0 ) )
[825]368
[5443]369                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
370                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
371               END DO
372            END DO
373         END DO
374      END DO
[825]375
[4688]376      tn_ice (:,:,:) = t_su (:,:,:)
377
378      ELSE 
[5443]379         ! if ln_iceini=false
[4688]380         a_i  (:,:,:) = 0._wp
381         v_i  (:,:,:) = 0._wp
382         v_s  (:,:,:) = 0._wp
383         smv_i(:,:,:) = 0._wp
384         oa_i (:,:,:) = 0._wp
385         ht_i (:,:,:) = 0._wp
386         ht_s (:,:,:) = 0._wp
387         sm_i (:,:,:) = 0._wp
388         o_i  (:,:,:) = 0._wp
389
390         e_i(:,:,:,:) = 0._wp
391         e_s(:,:,:,:) = 0._wp
392
393         DO jl = 1, jpl
394            DO jk = 1, nlay_i
[5443]395               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]396            END DO
397            DO jk = 1, nlay_s
[5443]398               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]399            END DO
400         END DO
401     
[5443]402      ENDIF ! ln_iceini
[4688]403     
404      at_i (:,:) = 0.0_wp
405      DO jl = 1, jpl
406         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
407      END DO
408      !
[825]409      !--------------------------------------------------------------------
[4161]410      ! 4) Global ice variables for output diagnostics                    |
[825]411      !--------------------------------------------------------------------
[4161]412      u_ice (:,:)     = 0._wp
413      v_ice (:,:)     = 0._wp
414      stress1_i(:,:)  = 0._wp
415      stress2_i(:,:)  = 0._wp
416      stress12_i(:,:) = 0._wp
[825]417
418      !--------------------------------------------------------------------
[4161]419      ! 5) Moments for advection
[825]420      !--------------------------------------------------------------------
421
[4161]422      sxopw (:,:) = 0._wp 
423      syopw (:,:) = 0._wp 
424      sxxopw(:,:) = 0._wp 
425      syyopw(:,:) = 0._wp 
426      sxyopw(:,:) = 0._wp
[3610]427
[4161]428      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
429      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
430      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
431      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
432      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
[825]433
[4161]434      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
435      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
436      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
437      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
438      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
[825]439
[4161]440      sxsal  (:,:,:)  = 0._wp
441      sysal  (:,:,:)  = 0._wp
442      sxxsal (:,:,:)  = 0._wp
443      syysal (:,:,:)  = 0._wp
444      sxysal (:,:,:)  = 0._wp
[825]445
[4335]446      sxage  (:,:,:)  = 0._wp
447      syage  (:,:,:)  = 0._wp
448      sxxage (:,:,:)  = 0._wp
449      syyage (:,:,:)  = 0._wp
450      sxyage (:,:,:)  = 0._wp
451
[825]452
[4688]453      CALL wrk_dealloc( jpi, jpj, zswitch )
[4161]454      CALL wrk_dealloc( jpi, jpj, zhemis )
[4688]455      CALL wrk_dealloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
456      CALL wrk_dealloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
[4161]457
[825]458   END SUBROUTINE lim_istate
459
460   SUBROUTINE lim_istate_init
461      !!-------------------------------------------------------------------
462      !!                   ***  ROUTINE lim_istate_init  ***
463      !!       
464      !! ** Purpose : Definition of initial state of the ice
465      !!
[4161]466      !! ** Method : Read the namiceini namelist and check the parameter
467      !!       values called at the first timestep (nit000)
[825]468      !!
[4161]469      !! ** input :
470      !!        Namelist namiceini
471      !!
472      !! history :
473      !!  8.5  ! 03-08 (C. Ethe) original code
474      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
[825]475      !!-----------------------------------------------------------------------------
[5443]476      NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s,  &
477         &                                      rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s
[4147]478      INTEGER :: ios                 ! Local integer output status for namelist read
[825]479      !!-----------------------------------------------------------------------------
[2528]480      !
[4147]481      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
482      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
483901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
484
485      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
486      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
487902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
[4624]488      IF(lwm) WRITE ( numoni, namiceini )
[4161]489
490      ! Define the initial parameters
491      ! -------------------------
492
493      IF(lwp) THEN
[825]494         WRITE(numout,*)
495         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
496         WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5443]497         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
498         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
499         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
500         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
501         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
502         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
503         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
504         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
505         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
506         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
507         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
508         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
[825]509      ENDIF
[4161]510
[825]511   END SUBROUTINE lim_istate_init
512
513#else
514   !!----------------------------------------------------------------------
515   !!   Default option :         Empty module          NO LIM sea-ice model
516   !!----------------------------------------------------------------------
517CONTAINS
518   SUBROUTINE lim_istate          ! Empty routine
519   END SUBROUTINE lim_istate
520#endif
521
522   !!======================================================================
523END MODULE limistate
Note: See TracBrowser for help on using the repository browser.