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 trunk/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90 @ 2864

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