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.
limsbc_2.F90 in branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90 @ 4148

Last change on this file since 4148 was 4148, checked in by cetlod, 10 years ago

merge in trunk changes between r3853 and r3940 and commit the changes, see ticket #1169

  • Property svn:keywords set to Id
File size: 27.9 KB
Line 
1MODULE limsbc_2
2   !!======================================================================
3   !!                       ***  MODULE limsbc_2   ***
4   !! LIM-2 :   updates the fluxes at the ocean surface with ice-ocean fluxes
5   !!======================================================================
6   !! History :  LIM  ! 2000-01 (H. Goosse) Original code
7   !!            1.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90
8   !!            3.0  ! 2006-07 (G. Madec) surface module
9   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
10   !!             -   ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step
11   !!           3.3.1 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation
12   !!            3.5  ! 2012-11 ((G. Madec, Y. Aksenov, A. Coward) salt and heat fluxes associated with e-p
13   !!----------------------------------------------------------------------
14#if defined key_lim2
15   !!----------------------------------------------------------------------
16   !!   'key_lim2'                                    LIM 2.0 sea-ice model
17   !!----------------------------------------------------------------------
18   !!   lim_sbc_alloc_2 : allocate the limsbc arrays
19   !!   lim_sbc_init    : initialisation
20   !!   lim_sbc_flx_2   : update mass, heat and salt fluxes at the ocean surface
21   !!   lim_sbc_tau_2   : update i- and j-stresses, and its modulus at the ocean surface
22   !!----------------------------------------------------------------------
23   USE par_oce          ! ocean parameters
24   USE phycst           ! physical constants
25   USE dom_oce          ! ocean domain
26   USE dom_ice_2        ! LIM-2: ice domain
27   USE ice_2            ! LIM-2: ice variables
28   USE sbc_ice          ! surface boundary condition: ice
29   USE sbc_oce          ! surface boundary condition: ocean
30   USE sbccpl
31   USE cpl_oasis3, ONLY : lk_cpl
32   USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
33   USE albedo           ! albedo parameters
34   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
35   USE lib_mpp          ! MPP library
36   USE wrk_nemo         ! work arrays
37   USE in_out_manager   ! I/O manager
38   USE diaar5, ONLY :   lk_diaar5
39   USE iom              ! I/O library
40   USE prtctl           ! Print control
41   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   lim_sbc_init_2     ! called by ice_init_2
47   PUBLIC   lim_sbc_flx_2      ! called by sbc_ice_lim_2
48   PUBLIC   lim_sbc_tau_2      ! called by sbc_ice_lim_2
49
50   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice
51   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
52   REAL(wp)  ::   rzero  = 0._wp       !     -      -
53   REAL(wp)  ::   rone   = 1._wp       !     -      -
54   !
55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0, sice_0   ! constant SSS and ice salinity used in levitating sea-ice case
56   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress              [N/m2]
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean relative velocity   [m/s]
58
59   !! * Substitutions
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/LIM2 4.0 , UCL - NEMO Consortium (2011)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   INTEGER FUNCTION lim_sbc_alloc_2()
69      !!-------------------------------------------------------------------
70      !!             ***  ROUTINE lim_sbc_alloc_2 ***
71      !!-------------------------------------------------------------------
72      ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) ,                       &
73         &      sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=lim_sbc_alloc_2)
74         !
75      IF( lk_mpp               )   CALL mpp_sum( lim_sbc_alloc_2 )
76      IF( lim_sbc_alloc_2 /= 0 )   CALL ctl_warn('lim_sbc_alloc_2: failed to allocate arrays.')
77      !
78   END FUNCTION lim_sbc_alloc_2
79
80
81   SUBROUTINE lim_sbc_flx_2( kt )
82      !!-------------------------------------------------------------------
83      !!                ***  ROUTINE lim_sbc_2 ***
84      !! 
85      !! ** Purpose :   Update surface ocean boundary condition over areas
86      !!              that are at least partially covered by sea-ice
87      !!         
88      !! ** Action  : - comput. of the momentum, heat and freshwater/salt
89      !!                fluxes at the ice-ocean interface.
90      !!              - Update the fluxes provided to the ocean
91      !!     
92      !! ** Outputs : - qsr     : sea heat flux    : solar
93      !!              - qns     : sea heat flux    : non solar (including heat content of the mass flux)
94      !!              - emp     : freshwater budget: mass flux
95      !!              - sfx     : freshwater budget: salt flux due to Freezing/Melting
96      !!              - utau    : sea surface i-stress (ocean referential)
97      !!              - vtau    : sea surface j-stress (ocean referential)
98      !!              - fr_i    : ice fraction
99      !!              - tn_ice  : sea-ice surface temperature
100      !!              - alb_ice : sea-ice alberdo (lk_cpl=T)
101      !!
102      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
103      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
104      !!---------------------------------------------------------------------
105      INTEGER, INTENT(in) ::   kt    ! number of iteration
106      !!
107      INTEGER  ::   ji, jj   ! dummy loop indices
108      INTEGER  ::   ii0, ii1, ij0, ij1         ! local integers
109      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
110      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
111      REAL(wp) ::   zqsr,     zqns,   zfmm     ! local scalars
112      REAL(wp) ::   zinda,    zfsalt, zemp     !   -      -
113      REAL(wp) ::   zemp_snw, zqhc,   zcd      !   -      -
114      REAL(wp) ::   zswitch                    !   -      -
115      REAL(wp), POINTER, DIMENSION(:,:)   ::   zqnsoce       ! 2D workspace
116      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace
117      !!---------------------------------------------------------------------
118     
119      CALL wrk_alloc( jpi, jpj, zqnsoce )
120      CALL wrk_alloc( jpi, jpj, 1, zalb, zalbp )
121
122      SELECT CASE( nn_ice_embd )                 ! levitating or embedded sea-ice option
123        CASE( 0    )   ;   zswitch = 1           ! (0) standard levitating sea-ice : salt exchange only
124        CASE( 1, 2 )   ;   zswitch = 0           ! (1) levitating sea-ice: salt and volume exchange but no pressure effect
125                                                 ! (2) embedded sea-ice : salt and volume fluxes and pressure
126      END SELECT                                 !   
127
128      !------------------------------------------!
129      !      heat flux at the ocean surface      !
130      !------------------------------------------!
131
132      zqnsoce(:,:) = qns(:,:)
133      DO jj = 1, jpj
134         DO ji = 1, jpi
135            zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
136            ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
137            i1mfr   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) )    ) )
138            idfr    = 1.0   - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
139            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
140            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
141            iadv    = ( 1  - i1mfr ) * zinda
142            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
143            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
144
145!!$            attempt to explain the tricky flags set above....
146!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo)
147!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice thermo)
148!!$
149!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = zinda if previous thermodynamic step overmelted the ice???
150!!$            ELSE                             ;   ifvt = 0.         !
151!!$            ENDIF
152!!$
153!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermodynamics
154!!$            ELSE                                     ;   idfr = 1.   
155!!$            ENDIF
156!!$
157!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous time and ice-free ocean currently
158!!$
159!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
160!!$                    = i1mfr if ifvt = 1 i.e.
161!!$                    = idfr  if ifvt = 0
162!!$!                 snow no ice   ice         ice or nothing  lead fraction increases
163!!$!                 at previous   now           at previous
164!!$!                -> ice area increases  ???         -> ice area decreases ???
165!!$
166!!$            iadv    = ( 1  - i1mfr ) * zinda
167!!$!                     pure ocean      ice at
168!!$!                     at current      previous
169!!$!                        -> = 1. if ice disapear between previous and current
170!!$
171!!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) ) 
172!!$!                            ice at     ???
173!!$!                            current         
174!!$!                         -> ???
175!!$
176!!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv
177!!$!                                                    ice disapear
178!!$
179!!$
180
181            !   computation the solar flux at ocean surface
182#if defined key_coupled 
183            zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
184#else
185            zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
186#endif           
187            !  computation the non solar heat flux at ocean surface
188            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr                                              &   ! part of the solar energy used in leads
189               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                             &
190               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice  &
191               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice 
192
193            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! store residual heat flux (to put into the ocean at the next time-step)
194            zqhc = ( rdq_snw(ji,jj)                                     &
195                 & + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice       ! heat flux due to snow ( & ice heat content,
196            !                                                            !           if ice/ocean mass exchange active)
197            qsr  (ji,jj) = zqsr                                          ! solar heat flux
198            qns  (ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux
199            !
200            !                          !------------------------------------------!
201            !                          !  mass and salt flux at the ocean surface !
202            !                          !------------------------------------------!
203            !
204            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
205#if defined key_coupled
206            !                                                  ! coupled mode:
207            zemp = + emp_tot(ji,jj)                            &     ! net mass flux over the grid cell (ice+ocean area)
208               &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )          ! minus the mass flux intercepted by sea-ice
209#else
210            !                                                  ! forced  mode:
211            zemp = + emp(ji,jj)     *         frld(ji,jj)      &     ! mass flux over open ocean fraction
212               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &     ! liquid precip. over ice reaches directly the ocean
213               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )          ! snow is intercepted by sea-ice (previous frld)
214#endif           
215            !
216            ! mass flux at the ocean/ice interface (sea ice fraction)
217            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                    ! snow melting = pure water that enters the ocean
218            zfmm     = rdm_ice(ji,jj) * r1_rdtice                    ! Freezing minus Melting (F-M)
219
220            fmmflx(ji,jj) = zfmm                                     ! F/M mass flux save at least for biogeochemical model
221
222            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
223            zfsalt = - sice_0(ji,jj) * zfmm                          ! F-M salt exchange
224            zcd    =   soce_0(ji,jj) * zfmm                          ! concentration/dilution term due to F-M
225            !
226            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux
227            ! salt and mass fluxes : non concentration dilution term in salt flux  and add F-M term in volume flux
228            sfx (ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange)
229            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (+ F/M mass flux if ice/ocean mass exchange)
230            !
231         END DO
232      END DO
233      !                                !------------------------------------------!
234      !                                !    mass of snow and ice per unit area    !
235      !                                !------------------------------------------!
236      IF( nn_ice_embd /= 0 ) THEN      ! embedded sea-ice (mass required)
237         snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step
238         !                                                      ! new mass per unit area
239         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) )
240         !                                                      ! time evolution of snow+ice mass
241         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / rdt_ice
242      ENDIF
243
244      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )     
245      CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )     
246      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
247
248      IF( lk_diaar5 ) THEN       ! AR5 diagnostics
249         CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * r1_rdtice )
250         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdm_ice(:,:) * r1_rdtice )
251         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice )
252      ENDIF
253
254      !-----------------------------------------------!
255      !   Coupling variables                          !
256      !-----------------------------------------------!
257
258#if defined key_coupled
259      tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
260      ht_i(:,:,1) = hicif(:,:)
261      ht_s(:,:,1) = hsnif(:,:)
262      a_i(:,:,1) = fr_i(:,:)
263      !                                  ! Computation of snow/ice and ocean albedo
264      CALL albedo_ice( tn_ice, ht_i, ht_s, zalbp, zalb )
265      alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
266      CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo
267#endif
268
269      IF(ln_ctl) THEN            ! control print
270         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
271         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=sfx   , clinfo2=' sfx     : ')
272         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
273            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
274         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
275      ENDIF 
276      !
277      CALL wrk_dealloc( jpi, jpj, zqnsoce )
278      CALL wrk_dealloc( jpi, jpj, 1, zalb, zalbp )
279      !
280   END SUBROUTINE lim_sbc_flx_2
281
282
283   SUBROUTINE lim_sbc_tau_2( kt , pu_oce, pv_oce )
284      !!-------------------------------------------------------------------
285      !!                ***  ROUTINE lim_sbc_tau_2 ***
286      !! 
287      !! ** Purpose : Update the ocean surface stresses due to the ice
288      !!         
289      !! ** Action  : * at each ice time step (every nn_fsbc time step):
290      !!                - compute the modulus of ice-ocean relative velocity
291      !!                  at T-point (C-grid) or I-point (B-grid)
292      !!                      tmod_io = rhoco * | U_ice-U_oce |
293      !!                - update the modulus of stress at ocean surface
294      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
295      !!              * at each ocean time step (each kt):
296      !!                  compute linearized ice-ocean stresses as
297      !!                      Utau = tmod_io * | U_ice - pU_oce |
298      !!                using instantaneous current ocean velocity (usually before)
299      !!
300      !!    NB: - the averaging operator used depends on the ice dynamics grid (cp_ice_msh='I' or 'C')
301      !!        - ice-ocean rotation angle only allowed in cp_ice_msh='I' case
302      !!        - here we make an approximation: taum is only computed every ice time step
303      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
304      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximation...
305      !!
306      !! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
307      !!              - taum       : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
308      !!---------------------------------------------------------------------
309      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
310      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
311      !!
312      INTEGER  ::   ji, jj   ! dummy loop indices
313      REAL(wp) ::   zfrldu, zat_u, zu_i, zutau_ice, zu_t, zmodt   ! local scalar
314      REAL(wp) ::   zfrldv, zat_v, zv_i, zvtau_ice, zv_t, zmodi   !   -      -
315      REAL(wp) ::   zsang, zumt                                   !    -         -
316      REAL(wp), POINTER, DIMENSION(:,:) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
317      !!---------------------------------------------------------------------
318      !
319      CALL wrk_alloc( jpi, jpj, ztio_u, ztio_v )
320      !
321      SELECT CASE( cp_ice_msh )     
322      !                             !-----------------------!
323      CASE( 'I' )                   !  B-grid ice dynamics  !   I-point (i.e. F-point with sea-ice indexation)
324         !                          !--=--------------------!
325         !
326         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
327!CDIR NOVERRCHK
328            DO jj = 1, jpj                               !* modulus of ice-ocean relative velocity at I-point
329!CDIR NOVERRCHK
330               DO ji = 1, jpi
331                  zu_i  = u_ice(ji,jj) - u_oce(ji,jj)                   ! ice-ocean relative velocity at I-point
332                  zv_i  = v_ice(ji,jj) - v_oce(ji,jj)
333                  tmod_io(ji,jj) = SQRT( zu_i * zu_i + zv_i * zv_i )    ! modulus of this velocity (at I-point)
334               END DO
335            END DO
336!CDIR NOVERRCHK
337            DO jj = 1, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
338!CDIR NOVERRCHK
339               DO ji = 1, jpim1   ! NO vector opt.
340                  !                                               ! modulus of U_ice-U_oce at T-point
341                  zumt  = 0.25_wp * (  tmod_io(ji+1,jj) + tmod_io(ji  ,jj  )    &   
342                     &               + tmod_io(ji,jj+1) + tmod_io(ji+1,jj+1)  )
343                  !                                               ! update the modulus of stress at ocean surface
344                  taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zumt * zumt
345               END DO
346            END DO
347            CALL lbc_lnk( taum, 'T', 1. )
348            !
349            utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
350            vtau_oce(:,:) = vtau(:,:)
351            !
352         ENDIF
353         !
354         !                                        !==  at each ocean time-step  ==!
355         !
356         !                                               !* ice/ocean stress WITH a ice-ocean rotation angle at I-point
357         DO jj = 2, jpj
358            zsang  = SIGN( 1._wp, gphif(1,jj) ) * sangvg          ! change the cosine angle sign in the SH
359            DO ji = 2, jpi    ! NO vect. opt. possible
360               ! ... ice-ocean relative velocity at I-point using instantaneous surface ocean current at u- & v-pts
361               zu_i = u_ice(ji,jj) - 0.5_wp * ( pu_oce(ji-1,jj  ) + pu_oce(ji-1,jj-1) ) * tmu(ji,jj)
362               zv_i = v_ice(ji,jj) - 0.5_wp * ( pv_oce(ji  ,jj-1) + pv_oce(ji-1,jj-1) ) * tmu(ji,jj)
363               ! ... components of stress with a ice-ocean rotation angle
364               zmodi = rhoco * tmod_io(ji,jj)                     
365               ztio_u(ji,jj) = zmodi * ( cangvg * zu_i - zsang * zv_i )
366               ztio_v(ji,jj) = zmodi * ( cangvg * zv_i + zsang * zu_i )
367            END DO
368         END DO
369         !                                               !* surface ocean stresses at u- and v-points
370         DO jj = 2, jpjm1
371            DO ji = 2, jpim1   ! NO vector opt.
372               !                                   ! ice-ocean stress at U and V-points  (from I-point values)
373               zutau_ice  = 0.5_wp * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
374               zvtau_ice  = 0.5_wp * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
375               !                                   ! open-ocean (lead) fraction at U- & V-points (from T-point values)
376               zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
377               zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
378               !                                   ! update the surface ocean stress (ice-cover wheighted)
379               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
380               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
381            END DO
382         END DO
383         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition
384         !
385         !
386         !                          !-----------------------!
387      CASE( 'C' )                   !  C-grid ice dynamics  !   U & V-points (same as in the ocean)
388         !                          !--=--------------------!
389         !
390         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
391!CDIR NOVERRCHK
392            DO jj = 2, jpjm1                          !* modulus of the ice-ocean velocity at T-point
393!CDIR NOVERRCHK
394               DO ji = fs_2, fs_jpim1
395                  zu_t  = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   ! 2*(U_ice-U_oce) at T-point
396                  zv_t  = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)     
397                  zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )                      ! |U_ice-U_oce|^2
398                  !                                               ! update the modulus of stress at ocean surface
399                  taum   (ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zmodt
400                  tmod_io(ji,jj) = SQRT( zmodt ) * rhoco          ! rhoco*|Uice-Uoce|
401               END DO
402            END DO
403            CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
404            !
405            utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step
406            vtau_oce(:,:) = vtau(:,:)
407            !
408         ENDIF
409         !
410         !                                        !==  at each ocean time-step  ==!
411         !
412         DO jj = 2, jpjm1                             !* ice stress over ocean WITHOUT a ice-ocean rotation angle
413            DO ji = fs_2, fs_jpim1
414               !                                            ! ocean area at u- & v-points
415               zfrldu  = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
416               zfrldv  = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
417               !                                            ! quadratic drag formulation without rotation
418               !                                            ! using instantaneous surface ocean current
419               zutau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
420               zvtau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
421               !                                            ! update the surface ocean stress (ice-cover wheighted)
422               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
423               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
424            END DO
425         END DO
426         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
427         !
428      END SELECT
429
430      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
431         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
432     
433      CALL wrk_dealloc( jpi, jpj, ztio_u, ztio_v )
434      !
435   END SUBROUTINE lim_sbc_tau_2
436
437
438   SUBROUTINE lim_sbc_init_2
439      !!-------------------------------------------------------------------
440      !!                  ***  ROUTINE lim_sbc_init  ***
441      !!             
442      !! ** Purpose : Preparation of the file ice_evolu for the output of
443      !!      the temporal evolution of key variables
444      !!
445      !! ** input   : Namelist namicedia
446      !!-------------------------------------------------------------------
447      !
448      IF(lwp) WRITE(numout,*)
449      IF(lwp) WRITE(numout,*) 'lim_sbc_init_2 : LIM-2 sea-ice - surface boundary condition'
450      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   '
451
452      !                                      ! allocate lim_sbc arrays
453      IF( lim_sbc_alloc_2() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_flx_2 : unable to allocate arrays' )
454      !
455      r1_rdtice = 1._wp / rdt_ice
456      !
457      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
458      sice_0(:,:) = sice
459      !
460      IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea
461         WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
462            &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
463            soce_0(:,:) = 4._wp
464            sice_0(:,:) = 2._wp
465         END WHERE
466      ENDIF
467      !                                      ! embedded sea ice
468      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
469         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) )
470         snwice_mass_b(:,:) = snwice_mass(:,:)
471      ELSE
472         snwice_mass  (:,:) = 0.e0           ! no mass exchanges
473         snwice_mass_b(:,:) = 0.e0           ! no mass exchanges
474      ENDIF
475      IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart :
476         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
477         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
478         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
479      ENDIF
480      !
481   END SUBROUTINE lim_sbc_init_2
482
483#else
484   !!----------------------------------------------------------------------
485   !!   Default option         Empty module        NO LIM 2.0 sea-ice model
486   !!----------------------------------------------------------------------
487#endif 
488
489   !!======================================================================
490END MODULE limsbc_2
Note: See TracBrowser for help on using the repository browser.