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 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

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