source: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limistate.F90 @ 2113

Last change on this file since 2113 was 2113, checked in by omamce, 11 years ago

O.M.

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