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/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • Property svn:keywords set to Id
File size: 25.5 KB
Line 
1MODULE limistate
2   !!======================================================================
3   !!                     ***  MODULE  limistate  ***
4   !!              Initialisation of diagnostics ice variables
5   !!======================================================================
6   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
7   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
8   !!----------------------------------------------------------------------
9#if defined key_lim3
10   !!----------------------------------------------------------------------
11   !!   'key_lim3' :                                    LIM3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   lim_istate      :  Initialisation of diagnostics ice variables
14   !!   lim_istate_init :  initialization of ice state and namelist read
15   !!----------------------------------------------------------------------
16   USE phycst           ! physical constant
17   USE oce              ! dynamics and tracers variables
18   USE dom_oce          ! ocean domain
19   USE sbc_oce          ! Surface boundary condition: ocean fields
20   USE eosbn2           ! equation of state
21   USE ice              ! sea-ice variables
22   USE par_ice          ! ice parameters
23   USE dom_ice          ! sea-ice domain
24   USE in_out_manager   ! I/O manager
25   USE lbclnk           ! lateral boundary condition - MPP exchanges
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_istate      ! routine called by lim_init.F90
31
32   !                                  !!** init namelist (namiceini) **
33   REAL(wp) ::   ttest    = 2.0_wp     ! threshold water temperature for initial sea ice
34   REAL(wp) ::   hninn    = 0.5_wp     ! initial snow thickness in the north
35   REAL(wp) ::   hginn_u  = 2.5_wp     ! initial ice thickness in the north
36   REAL(wp) ::   aginn_u  = 0.7_wp     ! initial leads area in the north
37   REAL(wp) ::   hginn_d  = 5.0_wp     ! initial ice thickness in the north
38   REAL(wp) ::   aginn_d  = 0.25_wp    ! initial leads area in the north
39   REAL(wp) ::   hnins    = 0.1_wp     ! initial snow thickness in the south
40   REAL(wp) ::   hgins_u  = 1.0_wp     ! initial ice thickness in the south
41   REAL(wp) ::   agins_u  = 0.7_wp     ! initial leads area in the south
42   REAL(wp) ::   hgins_d  = 2.0_wp     ! initial ice thickness in the south
43   REAL(wp) ::   agins_d  = 0.2_wp     ! initial leads area in the south
44   REAL(wp) ::   sinn     = 6.301_wp   ! initial salinity
45   REAL(wp) ::   sins     = 6.301_wp   !
46
47   !!----------------------------------------------------------------------
48   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE lim_istate
55      !!-------------------------------------------------------------------
56      !!                    ***  ROUTINE lim_istate  ***
57      !!
58      !! ** Purpose :   defined the sea-ice initial state
59      !!
60      !! ** Method  :   restart from a state defined in a binary file
61      !!                or from arbitrary sea-ice conditions
62      !!-------------------------------------------------------------------
63      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
64      USE wrk_nemo, ONLY:   wrk_1d_1, wrk_1d_2, wrk_1d_3, wrk_1d_4
65      USE wrk_nemo, ONLY:   zidto => wrk_2d_1   ! ice indicator
66      !
67      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
68      REAL(wp) ::   zeps6, zeps, ztmelts, epsi06   ! local scalars
69      REAL(wp) ::   zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs 
70      REAL(wp), POINTER, DIMENSION(:) ::   zgfactorn, zhin 
71      REAL(wp), POINTER, DIMENSION(:) ::   zgfactors, zhis
72      !--------------------------------------------------------------------
73
74      IF(  wrk_in_use(1, 1,2)  ) THEN
75         CALL ctl_stop( 'lim_istate : requested workspace arrays unavailable.' )   ;   RETURN
76      END IF
77      zgfactorn => wrk_1d_1(1:jpm)   ;   zhin => wrk_1d_3(1:jpm)   ! Set-up pointers to sub-arrays of workspaces
78      zgfactors => wrk_1d_2(1:jpm)   ;   zhis => wrk_1d_4(1:jpm)
79
80      !--------------------------------------------------------------------
81      ! 1) Preliminary things
82      !--------------------------------------------------------------------
83      epsi06 = 1.e-6_wp
84
85      CALL lim_istate_init     !  reading the initials parameters of the ice
86
87!!gm  in lim2  the initialisation if only done if required in the namelist :
88!!gm      IF( .NOT. ln_limini ) THEN
89!!gm  this should be added in lim3 namelist...
90
91      !--------------------------------------------------------------------
92      ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s)              |
93      !--------------------------------------------------------------------
94
95      IF(lwp) WRITE(numout,*)
96      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
97      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
98
99      t_bo(:,:) = tfreez( sn(:,:,1) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
100
101      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
102         DO ji = 1, jpi
103            IF( tn(ji,jj,1)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0.e0      ! no ice
104            ELSE                                             ;   zidto(ji,jj) = 1.e0      !    ice
105            ENDIF
106         END DO
107      END DO
108
109      t_bo(:,:) = t_bo(:,:) + rt0                          ! t_bo converted from Celsius to Kelvin (rt0 over land)
110
111      ! constants for heat contents
112      zeps   = 1.e-20_wp
113      zeps6  = 1.e-06_wp
114
115      ! zgfactor for initial ice distribution
116      zgfactorn(:) = 0._wp
117      zgfactors(:) = 0._wp
118
119      ! first ice type
120      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
121         zhin (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
122         zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp )
123         zhis (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
124         zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp )
125      END DO ! jl
126      zgfactorn(1) = aginn_u / zgfactorn(1)
127      zgfactors(1) = agins_u / zgfactors(1)
128
129      ! -------------
130      ! new distribution, polynom of second order, conserving area and volume
131      zh1 = 0._wp
132      zh2 = 0._wp
133      zh3 = 0._wp
134      DO jl = 1, jpl
135         zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
136         zh1 = zh1 + zh
137         zh2 = zh2 + zh * zh
138         zh3 = zh3 + zh * zh * zh
139      END DO
140      IF(lwp) WRITE(numout,*) ' zh1 : ', zh1
141      IF(lwp) WRITE(numout,*) ' zh2 : ', zh2
142      IF(lwp) WRITE(numout,*) ' zh3 : ', zh3
143
144      zvol = aginn_u * hginn_u
145      zare = aginn_u
146      IF( jpl >= 2 ) THEN
147         zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
148         zan = ( zare - zbn*zh1 ) / zh2
149      ENDIF
150
151      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
152      IF(lwp) WRITE(numout,*) ' zare: ', zare
153      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
154      IF(lwp) WRITE(numout,*) ' zan : ', zan 
155
156      zvol = agins_u * hgins_u
157      zare = agins_u
158      IF( jpl >= 2 ) THEN
159         zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
160         zas = ( zare - zbs*zh1 ) / zh2
161      ENDIF
162
163      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
164      IF(lwp) WRITE(numout,*) ' zare: ', zare
165      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
166      IF(lwp) WRITE(numout,*) ' zan : ', zan 
167
168      !end of new lines
169      ! -------------
170!!!
171      ! retour a LIMA_MEC
172      !     ! second ice type
173      !     zdummy  = hi_max(ice_cat_bounds(2,1)-1)
174      !     hi_max(ice_cat_bounds(2,1)-1) = 0.0
175
176      !     ! here to change !!!!
177      !     jm = 2
178      !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
179      !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
180      !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
181      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
182      !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0)
183      !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
184      !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
185      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
186      !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0)
187      !     END DO ! jl
188      !     zgfactorn(2) = aginn_d / zgfactorn(2)
189      !     zgfactors(2) = agins_d / zgfactors(2)
190      !     hi_max(ice_cat_bounds(2,1)-1) = zdummy
191      ! END retour a LIMA_MEC
192!!!
193
194!!gm  optimisation :  loop over the ice categories inside the ji, jj loop !!!
195
196      DO jj = 1, jpj
197         DO ji = 1, jpi
198
199            !--- Northern hemisphere
200            !----------------------------------------------------------------
201            IF( fcor(ji,jj) >= 0._wp ) THEN   
202
203               !-----------------------
204               ! Ice area / thickness
205               !-----------------------
206
207               IF ( jpl .EQ. 1) THEN ! one category
208
209                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
210                     a_i(ji,jj,jl)    = zidto(ji,jj) * aginn_u
211                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hginn_u
212                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
213                  END DO
214
215               ELSE ! several categories
216
217                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
218                     zhin(1)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
219                     a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* & 
220                        (zhin(1)-hginn_u)/2.0) , epsi06)
221                     ! new line
222                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) )
223                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(1) 
224                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
225                  END DO
226
227               ENDIF
228
229
230!!!
231               ! retour a LIMA_MEC
232               !              !ridged ice
233               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
234               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
235               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories
236               !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
237               !                 a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* &
238               !                                    (zhin(2)-hginn_d)/2.0) , epsi06)
239               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(2)
240               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
241               !              END DO
242               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
243
244               !              !rafted ice
245               !              jl = 6
246               !              a_i(ji,jj,jl)       = 0.0
247               !              ht_i(ji,jj,jl)      = 0.0
248               !              v_i(ji,jj,jl)       = 0.0
249               ! END retour a LIMA_MEC
250!!!
251
252               DO jl = 1, jpl
253
254                  !-------------
255                  ! Snow depth
256                  !-------------
257                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hninn
258                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
259
260                  !---------------
261                  ! Ice salinity
262                  !---------------
263                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sinn  + ( 1.0 - zidto(ji,jj) ) * 0.1
264                  smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
265
266                  !----------
267                  ! Ice age
268                  !----------
269                  o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
270                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
271
272                  !------------------------------
273                  ! Sea ice surface temperature
274                  !------------------------------
275
276                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
277
278                  !------------------------------------
279                  ! Snow temperature and heat content
280                  !------------------------------------
281
282                  DO jk = 1, nlay_s
283                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
284                     ! Snow energy of melting
285                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
286                     ! Change dimensions
287                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
288                     ! Multiply by volume, so that heat content in 10^9 Joules
289                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
290                        v_s(ji,jj,jl)  / nlay_s
291                  END DO !jk
292
293                  !-----------------------------------------------
294                  ! Ice salinities, temperature and heat content
295                  !-----------------------------------------------
296
297                  DO jk = 1, nlay_i
298                     t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
299                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1
300                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
301
302                     ! heat content per unit volume
303                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
304                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
305                        +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
306                        - rcp      * ( ztmelts - rtt ) &
307                        )
308
309                     ! Correct dimensions to avoid big values
310                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
311
312                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
313                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
314                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
315                        nlay_i
316                  END DO ! jk
317
318               END DO ! jl
319
320            ELSE ! on fcor
321
322               !--- Southern hemisphere
323               !----------------------------------------------------------------
324
325               !-----------------------
326               ! Ice area / thickness
327               !-----------------------
328
329               IF ( jpl .EQ. 1) THEN ! one category
330
331                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
332                     a_i(ji,jj,jl)    = zidto(ji,jj) * agins_u
333                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hgins_u
334                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
335                  END DO
336
337               ELSE ! several categories
338
339                  !level ice
340                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories
341
342                     zhis(1)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
343                     a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * & 
344                        (zhis(1)-hgins_u)/2.0) , epsi06 )
345                     ! new line square distribution volume conserving
346                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) )
347                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(1) 
348                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
349
350                  END DO ! jl
351
352               ENDIF
353
354!!!
355               ! retour a LIMA_MEC
356               !              !ridged ice
357               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
358               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
359               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories
360               !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
361               !                 a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2)   &
362               !                    &          * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 )
363               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(2)
364               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
365               !              END DO
366               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
367
368               !              !rafted ice
369               !              jl = 6
370               !              a_i(ji,jj,jl)       = 0.0
371               !              ht_i(ji,jj,jl)      = 0.0
372               !              v_i(ji,jj,jl)       = 0.0
373               ! END retour a LIMA_MEC
374!!!
375
376               DO jl = 1, jpl !over thickness categories
377
378                  !---------------
379                  ! Snow depth
380                  !---------------
381
382                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hnins
383                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
384
385                  !---------------
386                  ! Ice salinity
387                  !---------------
388
389                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sins  + ( 1.0 - zidto(ji,jj) ) * 0.1
390                  smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
391
392                  !----------
393                  ! Ice age
394                  !----------
395
396                  o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
397                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
398
399                  !------------------------------
400                  ! Sea ice surface temperature
401                  !------------------------------
402
403                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
404
405                  !----------------------------------
406                  ! Snow temperature / heat content
407                  !----------------------------------
408
409                  DO jk = 1, nlay_s
410                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
411                     ! Snow energy of melting
412                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
413                     ! Change dimensions
414                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
415                     ! Multiply by volume, so that heat content in 10^9 Joules
416                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
417                        v_s(ji,jj,jl)  / nlay_s
418                  END DO
419
420                  !---------------------------------------------
421                  ! Ice temperature, salinity and heat content
422                  !---------------------------------------------
423
424                  DO jk = 1, nlay_i
425                     t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
426                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1
427                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
428
429                     ! heat content per unit volume
430                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
431                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
432                        +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
433                        - rcp      * ( ztmelts - rtt ) &
434                        )
435
436                     ! Correct dimensions to avoid big values
437                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
438
439                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
440                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
441                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
442                        nlay_i
443                  END DO !jk
444
445               END DO ! jl
446
447            ENDIF ! on fcor
448
449         END DO
450      END DO
451
452      !--------------------------------------------------------------------
453      ! 3) Global ice variables for output diagnostics                    |
454      !--------------------------------------------------------------------
455
456      fsbbq (:,:)     = 0.e0
457      u_ice (:,:)     = 0.e0
458      v_ice (:,:)     = 0.e0
459      stress1_i(:,:)  = 0.0
460      stress2_i(:,:)  = 0.0
461      stress12_i(:,:) = 0.0
462
463      !--------------------------------------------------------------------
464      ! 4) Moments for advection
465      !--------------------------------------------------------------------
466
467      sxice (:,:,:)  = 0.e0   ;   sxsn (:,:,:)  = 0.e0   ;   sxa  (:,:,:)  = 0.e0
468      syice (:,:,:)  = 0.e0   ;   sysn (:,:,:)  = 0.e0   ;   sya  (:,:,:)  = 0.e0
469      sxxice(:,:,:)  = 0.e0   ;   sxxsn(:,:,:)  = 0.e0   ;   sxxa (:,:,:)  = 0.e0
470      syyice(:,:,:)  = 0.e0   ;   syysn(:,:,:)  = 0.e0   ;   syya (:,:,:)  = 0.e0
471      sxyice(:,:,:)  = 0.e0   ;   sxysn(:,:,:)  = 0.e0   ;   sxya (:,:,:)  = 0.e0
472
473      sxc0  (:,:,:)  = 0.e0   ;   sxe  (:,:,:,:)= 0.e0   
474      syc0  (:,:,:)  = 0.e0   ;   sye  (:,:,:,:)= 0.e0   
475      sxxc0 (:,:,:)  = 0.e0   ;   sxxe (:,:,:,:)= 0.e0   
476      syyc0 (:,:,:)  = 0.e0   ;   syye (:,:,:,:)= 0.e0   
477      sxyc0 (:,:,:)  = 0.e0   ;   sxye (:,:,:,:)= 0.e0   
478
479      sxsal  (:,:,:)  = 0.e0
480      sysal  (:,:,:)  = 0.e0
481      sxxsal (:,:,:)  = 0.e0
482      syysal (:,:,:)  = 0.e0
483      sxysal (:,:,:)  = 0.e0
484
485      !--------------------------------------------------------------------
486      ! 5) Lateral boundary conditions                                    |
487      !--------------------------------------------------------------------
488
489      DO jl = 1, jpl
490         CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. )
491         CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. )
492         CALL lbc_lnk( v_s(:,:,jl)  , 'T', 1. )
493         CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. )
494         CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. )
495         !
496         CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. )
497         CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. )
498         CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. )
499         CALL lbc_lnk( o_i(:,:,jl)  , 'T', 1. )
500         CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. )
501         DO jk = 1, nlay_s
502            CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. )
503            CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. )
504         END DO
505         DO jk = 1, nlay_i
506            CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. )
507            CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
508         END DO
509         !
510         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl)
511      END DO
512
513      CALL lbc_lnk( at_i , 'T', 1. )
514      at_i(:,:) = tms(:,:) * at_i(:,:)                       ! put 0 over land
515      !
516      CALL lbc_lnk( fsbbq  , 'T', 1. )
517      !
518      IF(  wrk_not_released(1, 1,2)  ) THEN
519         CALL ctl_stop( 'lim_istate : failed to release workspace arrays.' )
520      END IF
521      !
522   END SUBROUTINE lim_istate
523
524
525   SUBROUTINE lim_istate_init
526      !!-------------------------------------------------------------------
527      !!                   ***  ROUTINE lim_istate_init  ***
528      !!       
529      !! ** Purpose : Definition of initial state of the ice
530      !!
531      !! ** Method :   Read the namiceini namelist and check the parameter
532      !!             values called at the first timestep (nit000)
533      !!
534      !! ** input  :   namelist namiceini
535      !!-----------------------------------------------------------------------------
536      NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins,   &
537         &                hgins_u, agins_u, hgins_d, agins_d, sinn, sins
538      !!-----------------------------------------------------------------------------
539      !
540      REWIND ( numnam_ice )               ! Read Namelist namiceini
541      READ   ( numnam_ice , namiceini )
542      !
543      IF(lwp) THEN                        ! control print
544         WRITE(numout,*)
545         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
546         WRITE(numout,*) '~~~~~~~~~~~~~~~'
547         WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest
548         WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn
549         WRITE(numout,*) '   initial undef ice thickness in the north     hginn_u    = ', hginn_u
550         WRITE(numout,*) '   initial undef ice concentr. in the north     aginn_u    = ', aginn_u         
551         WRITE(numout,*) '   initial  def  ice thickness in the north     hginn_d    = ', hginn_d
552         WRITE(numout,*) '   initial  def  ice concentr. in the north     aginn_d    = ', aginn_d         
553         WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins 
554         WRITE(numout,*) '   initial undef ice thickness in the north     hgins_u    = ', hgins_u
555         WRITE(numout,*) '   initial undef ice concentr. in the north     agins_u    = ', agins_u         
556         WRITE(numout,*) '   initial  def  ice thickness in the north     hgins_d    = ', hgins_d
557         WRITE(numout,*) '   initial  def  ice concentr. in the north     agins_d    = ', agins_d         
558         WRITE(numout,*) '   initial  ice salinity       in the north     sinn       = ', sinn
559         WRITE(numout,*) '   initial  ice salinity       in the south     sins       = ', sins
560      ENDIF
561      !
562   END SUBROUTINE lim_istate_init
563
564#else
565   !!----------------------------------------------------------------------
566   !!   Default option :         Empty module          NO LIM sea-ice model
567   !!----------------------------------------------------------------------
568CONTAINS
569   SUBROUTINE lim_istate          ! Empty routine
570   END SUBROUTINE lim_istate
571#endif
572
573   !!======================================================================
574END MODULE limistate
Note: See TracBrowser for help on using the repository browser.