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

source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90 @ 5226

Last change on this file since 5226 was 5226, checked in by cetlod, 9 years ago

NEMOGCM_dev_r5204_CNRS_PISCES_dcy : minor bugs correction

  • Property svn:keywords set to Id
File size: 30.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 domvvl           ! ocean vertical scale factors
27   USE dom_ice_2        ! LIM-2: ice domain
28   USE ice_2            ! LIM-2: ice variables
29   USE sbc_ice          ! surface boundary condition: ice
30   USE sbc_oce          ! surface boundary condition: ocean
31   USE sbccpl
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 iom              ! I/O library
39   USE prtctl           ! Print control
40   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   lim_sbc_init_2     ! called by ice_init_2
46   PUBLIC   lim_sbc_flx_2      ! called by sbc_ice_lim_2
47   PUBLIC   lim_sbc_tau_2      ! called by sbc_ice_lim_2
48
49   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice
50   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
51   REAL(wp)  ::   rzero  = 0._wp       !     -      -
52   REAL(wp)  ::   rone   = 1._wp       !     -      -
53   !
54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0, sice_0   ! constant SSS and ice salinity used in levitating sea-ice case
55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress              [N/m2]
56   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean relative velocity   [m/s]
57
58   !! * Substitutions
59#  include "vectopt_loop_substitute.h90"
60#  include "domzgr_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      !!              - fr_i    : ice fraction
97      !!              - tn_ice  : sea-ice surface temperature
98      !!              - alb_ice : sea-ice albedo (lk_cpl=T)
99      !!
100      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
101      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
102      !!---------------------------------------------------------------------
103      INTEGER, INTENT(in) ::   kt    ! number of iteration
104      !!
105      INTEGER  ::   ji, jj   ! dummy loop indices
106      INTEGER  ::   ii0, ii1, ij0, ij1         ! local integers
107      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
108      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
109      REAL(wp) ::   zqsr,     zqns,   zfmm     ! local scalars
110      REAL(wp) ::   zinda,    zfsalt, zemp     !   -      -
111      REAL(wp) ::   zemp_snw, zqhc,   zcd      !   -      -
112      REAL(wp) ::   zswitch                    !   -      -
113      REAL(wp), POINTER, DIMENSION(:,:)   ::   zqnsoce       ! 2D workspace
114      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace
115      !!---------------------------------------------------------------------
116     
117      CALL wrk_alloc( jpi, jpj, zqnsoce )
118      CALL wrk_alloc( jpi, jpj, 1, zalb, zalbp )
119
120      SELECT CASE( nn_ice_embd )                 ! levitating or embedded sea-ice option
121        CASE( 0    )   ;   zswitch = 1           ! (0) standard levitating sea-ice : salt exchange only
122        CASE( 1, 2 )   ;   zswitch = 0           ! (1) levitating sea-ice: salt and volume exchange but no pressure effect
123                                                 ! (2) embedded sea-ice : salt and volume fluxes and pressure
124      END SELECT                                 !   
125
126      ! make calls for heat fluxes before it is modified
127      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) )   !     solar flux at ocean surface
128      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , qsr_ice(:,:,1) * pfrld(:,:) ) !     solar flux at ice surface
129      IF( l_trcdm2dc ) THEN
130        IF( iom_use('qsr_oce_mean') )   CALL iom_put( "qsr_oce_mean" , qsr_mean(:,:) * pfrld(:,:) )   !   daily mean solar flux at ocean surface
131        IF( iom_use('qsr_ice_mean') )   CALL iom_put( "qsr_ice_mean" , qsr_ice_mean(:,:,1) * pfrld(:,:) ) !  daily mean solar flux at ice surface
132      ENDIF
133      !------------------------------------------!
134      !      heat flux at the ocean surface      !
135      !------------------------------------------!
136
137      zqnsoce(:,:) = qns(:,:)
138      DO jj = 1, jpj
139         DO ji = 1, jpi
140            zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
141            ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
142            i1mfr   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) )    ) )
143            idfr    = 1.0   - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
144            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
145            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
146            iadv    = ( 1  - i1mfr ) * zinda
147            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
148            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
149
150!!$            attempt to explain the tricky flags set above....
151!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo)
152!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice thermo)
153!!$
154!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = zinda if previous thermodynamic step overmelted the ice???
155!!$            ELSE                             ;   ifvt = 0.         !
156!!$            ENDIF
157!!$
158!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermodynamics
159!!$            ELSE                                     ;   idfr = 1.   
160!!$            ENDIF
161!!$
162!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous time and ice-free ocean currently
163!!$
164!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
165!!$                    = i1mfr if ifvt = 1 i.e.
166!!$                    = idfr  if ifvt = 0
167!!$!                 snow no ice   ice         ice or nothing  lead fraction increases
168!!$!                 at previous   now           at previous
169!!$!                -> ice area increases  ???         -> ice area decreases ???
170!!$
171!!$            iadv    = ( 1  - i1mfr ) * zinda
172!!$!                     pure ocean      ice at
173!!$!                     at current      previous
174!!$!                        -> = 1. if ice disapear between previous and current
175!!$
176!!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) ) 
177!!$!                            ice at     ???
178!!$!                            current         
179!!$!                         -> ???
180!!$
181!!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv
182!!$!                                                    ice disapear
183!!$
184!!$
185
186            !   computation the solar flux at ocean surface
187            IF( lk_cpl ) THEN
188               zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
189            ELSE
190               zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
191            ENDIF
192            !  computation the non solar heat flux at ocean surface
193            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr                                              &   ! part of the solar energy used in leads
194               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                             &
195               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice  &
196               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice 
197
198            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! store residual heat flux (to put into the ocean at the next time-step)
199            zqhc = ( rdq_snw(ji,jj)                                     &
200                 & + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice       ! heat flux due to snow ( & ice heat content,
201            !                                                            !           if ice/ocean mass exchange active)
202            qsr  (ji,jj) = zqsr                                          ! solar heat flux
203            qns  (ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux
204            !
205            !                          !------------------------------------------!
206            !                          !  mass and salt flux at the ocean surface !
207            !                          !------------------------------------------!
208            !
209            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
210            !                                                  ! coupled mode:
211            IF( lk_cpl ) THEN
212               zemp = + emp_tot(ji,jj)                            &     ! net mass flux over the grid cell (ice+ocean area)
213                  &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )          ! minus the mass flux intercepted by sea-ice
214            ELSE
215               !                                                  ! forced  mode:
216               zemp = + emp(ji,jj)     *         frld(ji,jj)      &     ! mass flux over open ocean fraction
217                  &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &     ! liquid precip. over ice reaches directly the ocean
218                  &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )          ! snow is intercepted by sea-ice (previous frld)
219            ENDIF
220            !
221            ! mass flux at the ocean/ice interface (sea ice fraction)
222            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                    ! snow melting = pure water that enters the ocean
223            zfmm     = rdm_ice(ji,jj) * r1_rdtice                    ! Freezing minus Melting (F-M)
224
225            fmmflx(ji,jj) = zfmm                                     ! F/M mass flux save at least for biogeochemical model
226
227            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
228            zfsalt = - sice_0(ji,jj) * zfmm                          ! F-M salt exchange
229            zcd    =   soce_0(ji,jj) * zfmm                          ! concentration/dilution term due to F-M
230            !
231            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux
232            ! salt and mass fluxes : non concentration dilution term in salt flux  and add F-M term in volume flux
233            sfx (ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange)
234            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (+ F/M mass flux if ice/ocean mass exchange)
235            !
236         END DO
237      END DO
238      !                                !------------------------------------------!
239      !                                !    mass of snow and ice per unit area    !
240      !                                !------------------------------------------!
241      IF( nn_ice_embd /= 0 ) THEN      ! embedded sea-ice (mass required)
242         snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step
243         !                                                      ! new mass per unit area
244         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) )
245         !                                                      ! time evolution of snow+ice mass
246         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / rdt_ice
247      ENDIF
248
249      IF( iom_use('hflx_ice_cea' ) )   CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )     
250      IF( iom_use('qns_io_cea'   ) )   CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )     
251      IF( iom_use('qsr_io_cea'   ) )   CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
252
253      IF( iom_use('isnwmlt_cea'  ) )   CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * r1_rdtice )
254      IF( iom_use('fsal_virt_cea') )   CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdm_ice(:,:) * r1_rdtice )
255      IF( iom_use('fsal_real_cea') )   CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice )
256
257      !-----------------------------------------------!
258      !   Coupling variables                          !
259      !-----------------------------------------------!
260
261      IF( lk_cpl) THEN
262         tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
263         ht_i(:,:,1) = hicif(:,:)
264         ht_s(:,:,1) = hsnif(:,:)
265         a_i(:,:,1) = fr_i(:,:)
266         !                                  ! Computation of snow/ice and ocean albedo
267         CALL albedo_ice( tn_ice, ht_i, ht_s, zalbp, zalb )
268         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
269         IF( iom_use('icealb_cea' ) )   CALL iom_put( 'icealb_cea', alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo
270      ENDIF
271
272      !   daily mean qsr when diurnal cycle is applied on physics - for BGC models
273      IF( l_trcdm2dc ) THEN 
274         !   computation the solar flux at ocean surface
275         IF( lk_cpl ) THEN
276            qsr_mean(:,:) = qsr_mean(:,:) + ( fstric_mean(:,:) - qsr_ice_mean(:,:,1) ) * ( 1.0 - pfrld(:,:) )  ! qsr_mean = qsr_tot
277         ELSE
278            qsr_mean(:,:) =  pfrld(:,:) * qsr_mean(:,:) + ( 1. - pfrld(:,:) ) * fstric_mean(:,:)
279         ENDIF
280      ENDIF
281
282
283      IF(ln_ctl) THEN            ! control print
284         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
285         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=sfx   , clinfo2=' sfx     : ')
286         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
287      ENDIF 
288      !
289      CALL wrk_dealloc( jpi, jpj, zqnsoce )
290      CALL wrk_dealloc( jpi, jpj, 1, zalb, zalbp )
291      !
292   END SUBROUTINE lim_sbc_flx_2
293
294
295   SUBROUTINE lim_sbc_tau_2( kt , pu_oce, pv_oce )
296      !!-------------------------------------------------------------------
297      !!                ***  ROUTINE lim_sbc_tau_2 ***
298      !! 
299      !! ** Purpose : Update the ocean surface stresses due to the ice
300      !!         
301      !! ** Action  : * at each ice time step (every nn_fsbc time step):
302      !!                - compute the modulus of ice-ocean relative velocity
303      !!                  at T-point (C-grid) or I-point (B-grid)
304      !!                      tmod_io = rhoco * | U_ice-U_oce |
305      !!                - update the modulus of stress at ocean surface
306      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
307      !!              * at each ocean time step (each kt):
308      !!                  compute linearized ice-ocean stresses as
309      !!                      Utau = tmod_io * | U_ice - pU_oce |
310      !!                using instantaneous current ocean velocity (usually before)
311      !!
312      !!    NB: - the averaging operator used depends on the ice dynamics grid (cp_ice_msh='I' or 'C')
313      !!        - ice-ocean rotation angle only allowed in cp_ice_msh='I' case
314      !!        - here we make an approximation: taum is only computed every ice time step
315      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
316      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximation...
317      !!
318      !! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
319      !!              - taum       : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
320      !!---------------------------------------------------------------------
321      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
322      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
323      !!
324      INTEGER  ::   ji, jj   ! dummy loop indices
325      REAL(wp) ::   zfrldu, zat_u, zu_i, zutau_ice, zu_t, zmodt   ! local scalar
326      REAL(wp) ::   zfrldv, zat_v, zv_i, zvtau_ice, zv_t, zmodi   !   -      -
327      REAL(wp) ::   zsang, zumt                                   !    -         -
328      REAL(wp), POINTER, DIMENSION(:,:) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
329      !!---------------------------------------------------------------------
330      !
331      CALL wrk_alloc( jpi, jpj, ztio_u, ztio_v )
332      !
333      SELECT CASE( cp_ice_msh )     
334      !                             !-----------------------!
335      CASE( 'I' )                   !  B-grid ice dynamics  !   I-point (i.e. F-point with sea-ice indexation)
336         !                          !--=--------------------!
337         !
338         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
339!CDIR NOVERRCHK
340            DO jj = 1, jpj                               !* modulus of ice-ocean relative velocity at I-point
341!CDIR NOVERRCHK
342               DO ji = 1, jpi
343                  zu_i  = u_ice(ji,jj) - u_oce(ji,jj)                   ! ice-ocean relative velocity at I-point
344                  zv_i  = v_ice(ji,jj) - v_oce(ji,jj)
345                  tmod_io(ji,jj) = SQRT( zu_i * zu_i + zv_i * zv_i )    ! modulus of this velocity (at I-point)
346               END DO
347            END DO
348!CDIR NOVERRCHK
349            DO jj = 1, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
350!CDIR NOVERRCHK
351               DO ji = 1, jpim1   ! NO vector opt.
352                  !                                               ! modulus of U_ice-U_oce at T-point
353                  zumt  = 0.25_wp * (  tmod_io(ji+1,jj) + tmod_io(ji  ,jj  )    &   
354                     &               + tmod_io(ji,jj+1) + tmod_io(ji+1,jj+1)  )
355                  !                                               ! update the modulus of stress at ocean surface
356                  taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zumt * zumt
357               END DO
358            END DO
359            CALL lbc_lnk( taum, 'T', 1. )
360            !
361            utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
362            vtau_oce(:,:) = vtau(:,:)
363            !
364         ENDIF
365         !
366         !                                        !==  at each ocean time-step  ==!
367         !
368         !                                               !* ice/ocean stress WITH a ice-ocean rotation angle at I-point
369         DO jj = 2, jpj
370            zsang  = SIGN( 1._wp, gphif(1,jj) ) * sangvg          ! change the cosine angle sign in the SH
371            DO ji = 2, jpi    ! NO vect. opt. possible
372               ! ... ice-ocean relative velocity at I-point using instantaneous surface ocean current at u- & v-pts
373               zu_i = u_ice(ji,jj) - 0.5_wp * ( pu_oce(ji-1,jj  ) + pu_oce(ji-1,jj-1) ) * tmu(ji,jj)
374               zv_i = v_ice(ji,jj) - 0.5_wp * ( pv_oce(ji  ,jj-1) + pv_oce(ji-1,jj-1) ) * tmu(ji,jj)
375               ! ... components of stress with a ice-ocean rotation angle
376               zmodi = rhoco * tmod_io(ji,jj)                     
377               ztio_u(ji,jj) = zmodi * ( cangvg * zu_i - zsang * zv_i )
378               ztio_v(ji,jj) = zmodi * ( cangvg * zv_i + zsang * zu_i )
379            END DO
380         END DO
381         !                                               !* surface ocean stresses at u- and v-points
382         DO jj = 2, jpjm1
383            DO ji = 2, jpim1   ! NO vector opt.
384               !                                   ! ice-ocean stress at U and V-points  (from I-point values)
385               zutau_ice  = 0.5_wp * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
386               zvtau_ice  = 0.5_wp * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
387               !                                   ! open-ocean (lead) fraction at U- & V-points (from T-point values)
388               zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
389               zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
390               !                                   ! update the surface ocean stress (ice-cover wheighted)
391               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
392               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
393            END DO
394         END DO
395         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition
396         !
397         !
398         !                          !-----------------------!
399      CASE( 'C' )                   !  C-grid ice dynamics  !   U & V-points (same as in the ocean)
400         !                          !--=--------------------!
401         !
402         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
403!CDIR NOVERRCHK
404            DO jj = 2, jpjm1                          !* modulus of the ice-ocean velocity at T-point
405!CDIR NOVERRCHK
406               DO ji = fs_2, fs_jpim1
407                  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
408                  zv_t  = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)     
409                  zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )                      ! |U_ice-U_oce|^2
410                  !                                               ! update the modulus of stress at ocean surface
411                  taum   (ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zmodt
412                  tmod_io(ji,jj) = SQRT( zmodt ) * rhoco          ! rhoco*|Uice-Uoce|
413               END DO
414            END DO
415            CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
416            !
417            utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step
418            vtau_oce(:,:) = vtau(:,:)
419            !
420         ENDIF
421         !
422         !                                        !==  at each ocean time-step  ==!
423         !
424         DO jj = 2, jpjm1                             !* ice stress over ocean WITHOUT a ice-ocean rotation angle
425            DO ji = fs_2, fs_jpim1
426               !                                            ! ocean area at u- & v-points
427               zfrldu  = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
428               zfrldv  = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
429               !                                            ! quadratic drag formulation without rotation
430               !                                            ! using instantaneous surface ocean current
431               zutau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
432               zvtau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
433               !                                            ! update the surface ocean stress (ice-cover wheighted)
434               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
435               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
436            END DO
437         END DO
438         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
439         !
440      END SELECT
441
442      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
443         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
444     
445      CALL wrk_dealloc( jpi, jpj, ztio_u, ztio_v )
446      !
447   END SUBROUTINE lim_sbc_tau_2
448
449
450   SUBROUTINE lim_sbc_init_2
451      !!-------------------------------------------------------------------
452      !!                  ***  ROUTINE lim_sbc_init  ***
453      !!             
454      !! ** Purpose : Preparation of the file ice_evolu for the output of
455      !!      the temporal evolution of key variables
456      !!
457      !! ** input   : Namelist namicedia
458      !!-------------------------------------------------------------------
459      !
460      INTEGER :: jk           ! local integer
461      !
462      IF(lwp) WRITE(numout,*)
463      IF(lwp) WRITE(numout,*) 'lim_sbc_init_2 : LIM-2 sea-ice - surface boundary condition'
464      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   '
465
466      !                                      ! allocate lim_sbc arrays
467      IF( lim_sbc_alloc_2() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_flx_2 : unable to allocate arrays' )
468      !
469      r1_rdtice = 1._wp / rdt_ice
470      !
471      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
472      sice_0(:,:) = sice
473      !
474      IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea
475         WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
476            &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
477            soce_0(:,:) = 4._wp
478            sice_0(:,:) = 2._wp
479         END WHERE
480      ENDIF
481      !                                      ! embedded sea ice
482      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
483         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) )
484         snwice_mass_b(:,:) = snwice_mass(:,:)
485      ELSE
486         snwice_mass  (:,:) = 0.e0           ! no mass exchanges
487         snwice_mass_b(:,:) = 0.e0           ! no mass exchanges
488         snwice_fmass (:,:) = 0.e0           ! no mass exchanges
489      ENDIF
490      IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart :
491         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
492         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
493         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
494         do jk = 1,jpkm1                     ! adjust initial vertical scale factors
495          fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
496          fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
497         end do
498         fse3t_a(:,:,:) = fse3t_b(:,:,:)
499         ! Reconstruction of all vertical scale factors at now and before time steps
500         ! =============================================================================
501         ! Horizontal scale factor interpolations
502         ! --------------------------------------
503         CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
504         CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
505         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
506         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
507         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
508         ! Vertical scale factor interpolations
509         ! ------------------------------------
510         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
511         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
512         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
513         CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
514         CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
515         ! t- and w- points depth
516         ! ----------------------
517         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
518         fsdepw_n(:,:,1) = 0.0_wp
519         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
520         DO jk = 2, jpk
521            fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
522            fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
523            fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
524         END DO
525      ENDIF
526      !
527   END SUBROUTINE lim_sbc_init_2
528
529#else
530   !!----------------------------------------------------------------------
531   !!   Default option         Empty module        NO LIM 2.0 sea-ice model
532   !!----------------------------------------------------------------------
533#endif 
534
535   !!======================================================================
536END MODULE limsbc_2
Note: See TracBrowser for help on using the repository browser.