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.F90 in branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90 @ 3977

Last change on this file since 3977 was 3977, checked in by flavoni, 11 years ago

add print of fwb correction value, and hminrhg for LIM2, in CNRS LIM3 branch

  • Property svn:keywords set to Id
File size: 25.3 KB
Line 
1MODULE limsbc
2   !!======================================================================
3   !!                       ***  MODULE limsbc   ***
4   !!           computation of the flux at the sea ice/ocean interface
5   !!======================================================================
6   !! History :   -   ! 2006-07 (M. Vancoppelle)  LIM3 original code
7   !!            3.0  ! 2008-03 (C. Tallandier)  surface module
8   !!             -   ! 2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling
9   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
10   !!                 !                  + simplification of the ice-ocean stress calculation
11   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
12   !!             -   ! 2012    (D. Iovino) salt flux change
13   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
14   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass
15   !!----------------------------------------------------------------------
16#if defined key_lim3
17   !!----------------------------------------------------------------------
18   !!   'key_lim3'                                    LIM 3.0 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   lim_sbc_alloc : allocate the limsbc arrays
21   !!   lim_sbc_init  : initialisation
22   !!   lim_sbc_flx   : updates mass, heat and salt fluxes at the ocean surface
23   !!   lim_sbc_tau   : update i- and j-stresses, and its modulus at the ocean surface
24   !!----------------------------------------------------------------------
25   USE par_oce          ! ocean parameters
26   USE par_ice          ! ice parameters
27   USE dom_oce          ! ocean domain
28   USE sbc_ice          ! Surface boundary condition: sea-ice fields
29   USE sbc_oce          ! Surface boundary condition: ocean fields
30   USE phycst           ! physical constants
31   USE albedo           ! albedo parameters
32   USE ice              ! LIM sea-ice variables
33   USE lbclnk           ! ocean lateral boundary condition
34   USE in_out_manager   ! I/O manager
35   USE lib_mpp          ! MPP library
36   USE wrk_nemo         ! work arrays
37   USE prtctl           ! Print control
38   USE cpl_oasis3, ONLY : lk_cpl
39   USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget
40   USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
41   USE dom_ice,    ONLY : tms
42   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   lim_sbc_init   ! called by ice_init
48   PUBLIC   lim_sbc_flx    ! called by sbc_ice_lim
49   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim
50
51   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice
52   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
53   REAL(wp)  ::   epsi20 = 1.e-20_wp   ! constant values
54   REAL(wp)  ::   rzero  = 0._wp   
55   REAL(wp)  ::   rone   = 1._wp
56   REAL(wp) ::   r1_rau0   ! = 1 / rau0
57
58   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
59   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
60   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0  , sice_0     ! cst SSS and ice salinity (levitating sea-ice)
61
62   !! * Substitutions
63#  include "vectopt_loop_substitute.h90"
64   !!----------------------------------------------------------------------
65   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
66   !! $Id$
67   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69CONTAINS
70
71   INTEGER FUNCTION lim_sbc_alloc()
72      !!-------------------------------------------------------------------
73      !!             ***  ROUTINE lim_sbc_alloc ***
74      !!-------------------------------------------------------------------
75      ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) ,                       &
76         &      sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=lim_sbc_alloc)
77         !
78      IF( lk_mpp             )   CALL mpp_sum( lim_sbc_alloc )
79      IF( lim_sbc_alloc /= 0 )   CALL ctl_warn('lim_sbc_alloc: failed to allocate arrays')
80   END FUNCTION lim_sbc_alloc
81
82
83   SUBROUTINE lim_sbc_flx( kt )
84      !!-------------------------------------------------------------------
85      !!                ***  ROUTINE lim_sbc_flx ***
86      !! 
87      !! ** Purpose :   Update the surface ocean boundary condition for heat
88      !!              salt and mass over areas where sea-ice is non-zero
89      !!         
90      !! ** Action  : - computes the heat and freshwater/salt fluxes
91      !!              at the ice-ocean interface.
92      !!              - Update the ocean sbc
93      !!     
94      !! ** Outputs : - qsr     : sea heat flux:     solar
95      !!              - qns     : sea heat flux: non solar
96      !!              - emp     : freshwater budget: volume flux
97      !!              - emps    : freshwater budget: concentration/dillution
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  ::   ierr             ! local integer
109      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
110      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
111      REAL(wp) ::   zinda, zindb, zfons, zpme              ! local scalars
112      REAL(wp) ::   zfmm             ! IOVINO freezing minus melting (F-M)
113      REAL(wp), POINTER, DIMENSION(:,:) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes
114      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace
115      !!---------------------------------------------------------------------
116     
117      CALL wrk_alloc( jpi, jpj, zfcm1 , zfcm2 )
118      IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
119
120      !------------------------------------------!
121      !      heat flux at the ocean surface      !
122      !------------------------------------------!
123      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
124      ! changed to old_frld and old ht_i
125
126      DO jj = 1, jpj
127         DO ji = 1, jpi
128            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
129            zindb   = 1.0 - MAX( rzero , SIGN( rone , - iatte(ji,jj) ) )
130            ifvt    = zinda  *  MAX( rzero , SIGN( rone, - phicif(ji,jj) ) )  !subscripts are bad here
131            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - at_i(ji,jj) ) )
132            idfr    = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(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 
138
139            ! switch --- 1.0 ---------------- 0.0 --------------------
140            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
141            ! zinda   | if pfrld = 1       | if pfrld < 1            |
142            !  -> ifvt| if pfrld old_ht_i
143            ! i1mfr   | if frld = 1        | if frld  < 1            |
144            ! idfr    | if frld <= pfrld    | if frld > pfrld        |
145            ! iflt    |
146            ! ial     |
147            ! iadv    |
148            ! ifral
149            ! ifrdv
150
151            !   computation the solar flux at ocean surface
152            zfcm1(ji,jj)   = pfrld(ji,jj) * qsr(ji,jj)  + &
153                 &           zindb * ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) / MAX( iatte(ji,jj), epsi20 )
154            ! fstric     Solar flux transmitted trough the ice
155            ! qsr        Net short wave heat flux on free ocean
156            ! new line
157            fscmbq(ji,jj) = zindb * ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) / MAX( iatte(ji,jj), epsi20 )
158
159            !  computation the non solar heat flux at ocean surface
160            zfcm2(ji,jj) = - zfcm1(ji,jj)                  &
161               &           + iflt    * ( fscmbq(ji,jj) )   & ! total abl -> fscmbq is given to the ocean
162               ! fscmbq and ffltbif are obsolete
163               !              &           + iflt * ffltbif(ji,jj) !!! only if one category is used
164               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   &
165               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice   &
166               &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!!
167               &           + fheat_rpo(ji,jj) & ! contribution from ridge formation
168               &           + fheat_res(ji,jj)
169            ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean
170            !         computed in limthd_zdf.F90
171            ! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t
172            ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok)
173            ! qldif   heat balance of the lead (or of the open ocean)
174            ! qfvbq   i think this is wrong!
175            ! ---> Array used to store energy in case of total lateral ablation
176            ! qfvbq latent heat uptake/release after accretion/ablation
177            ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
178
179            IF ( num_sal == 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + fhbri(ji,jj) ! new contribution due to brine drainage
180
181            ! bottom radiative component is sent to the computation of the
182            ! oceanic heat flux
183            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     
184
185            ! used to compute the oceanic heat flux at the next time step
186            qsr(ji,jj) = zfcm1(ji,jj)                                       ! solar heat flux
187            qns(ji,jj) = zfcm2(ji,jj) - fdtcn(ji,jj)                        ! non solar heat flux
188            !                           ! fdtcn : turbulent oceanic heat flux
189
190            !!gm   this IF prevents the vertorisation of the whole loop
191          !  IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN
192          !     WRITE(numout,*) 'lim_sbc : heat fluxes '
193          !     WRITE(numout,*) ' at_i      : ', at_i(jiindx,jjindx)
194          !     WRITE(numout,*) ' ht_i      : ', SUM( ht_i(jiindx,jjindx,1:jpl) )
195          !     WRITE(numout,*) ' ht_s      : ', SUM( ht_s(jiindx,jjindx,1:jpl) )
196          !     WRITE(numout,*)
197          !     WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx)
198          !     WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
199          !     WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx)
200          !     WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx)
201          !     WRITE(numout,*)
202          !     WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx)
203          !     WRITE(numout,*) ' zfcm2     : ', zfcm2(jiindx,jjindx)
204          !     WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
205          !     WRITE(numout,*) ' ifral     : ', ifral
206          !     WRITE(numout,*) ' ial       : ', ial 
207          !     WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx)
208          !     WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx)
209          !     !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice
210          !     !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice
211          !     WRITE(numout,*) ' ifrdv     : ', ifrdv
212          !     WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx)
213          !     WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx)
214          !     !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice
215          !     !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice
216          !     WRITE(numout,*) ' '
217          !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
218          !     WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx)
219          !     WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx)
220          !     WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx)
221          !     WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
222          !  ENDIF
223            !!gm   end
224         END DO
225      END DO
226
227      !------------------------------------------!
228      !      mass flux at the ocean surface      !
229      !------------------------------------------!
230
231!!gm   optimisation: this loop have to be merged with the previous one
232      DO jj = 1, jpj
233         DO ji = 1, jpi
234            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
235            !  -------------------------------------------------------------------------------------
236            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
237            !  Thus  FW  flux  =  External ( E-P+snow melt)
238            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
239            !                     Associated to Ice formation AND Ice melting
240            !                     Even if i see Ice melting as a FW and SALT flux
241            !       
242
243            !  computing freshwater exchanges at the ice/ocean interface
244            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
245               &   + tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
246               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
247               &   - rdmsnif(ji,jj) * r1_rdtice                       &   ! freshwaterflux due to snow melting
248               &   + fmmec(ji,jj)                                         ! snow falling when ridging
249
250
251            !  computing salt exchanges at the ice/ocean interface
252            !  sice should be the same as computed with the ice model
253            !zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice
254            ! SOCE
255            !zfons =  ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice
256            zfmm = rdmicif(ji,jj) * r1_rdtice  ! IOVINO 
257            !CT useless            !  salt flux for constant salinity
258            !CT useless            fsalt(ji,jj)      =  zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj)
259            !  salt flux for variable salinity
260            zinda             = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
261            !  correcting brine and salt fluxes
262            fsbri(ji,jj)      =  zinda*fsbri(ji,jj)
263            !  converting the salt fluxes from ice to a freshwater flux from ocean
264            ! fsalt_res(ji,jj)  =  fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 )
265            ! fseqv(ji,jj)      =  fseqv(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
266            ! fsbri(ji,jj)      =  fsbri(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
267            ! fsalt_rpo(ji,jj)  =  fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 )
268
269            !  freshwater mass exchange (positive to the ice, negative for the ocean ?)
270            !  actually it's a salt flux (so it's minus freshwater flux)
271            !  if sea ice grows, zfons is positive, fsalt also
272            !  POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN
273            !  POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1]
274
275            emp(ji,jj) =  - zpme + zfmm ! volume flux IOVINO 
276            ! emp(ji,jj) = - zpme
277         END DO
278      END DO
279
280      IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux
281         emps(:,:) =   fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ! + emp(:,:) ! IOVINO
282      ELSE                         ! constant ice salinity:
283         emps(:,:) =   fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ! + emp(:,:)              ! IOVINO
284      ENDIF
285
286      !-----------------------------------------------!
287      !   mass of snow and ice per unit area          !
288      !-----------------------------------------------!
289     !SF IF( nn_ice_embd /= 0 ) THEN                               ! embedded sea-ice (mass required)
290         snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step
291         !                                                      ! new mass per unit area
292         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
293         !                                                      ! time evolution of snow+ice mass
294         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
295      !SF ENDIF
296
297      !-----------------------------------------------!
298      !   Storing the transmitted variables           !
299      !-----------------------------------------------!
300      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
301      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
302
303      !------------------------------------------------!
304      !    Computation of snow/ice and ocean albedo    !
305      !------------------------------------------------!
306      IF( lk_cpl ) THEN          ! coupled case
307         CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )                  ! snow/ice albedo
308         !
309         alb_ice(:,:,:) =  0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys)
310      ENDIF
311
312      IF(ln_ctl) THEN
313         CALL prt_ctl( tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns , clinfo2=' qns     : ' )
314         CALL prt_ctl( tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps, clinfo2=' emps    : ' )
315         CALL prt_ctl( tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ' )
316         CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
317      ENDIF
318      !
319      CALL wrk_dealloc( jpi, jpj, zfcm1 , zfcm2 )
320      IF( lk_cpl )   CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp )
321      !
322   END SUBROUTINE lim_sbc_flx
323
324
325   SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce )
326      !!-------------------------------------------------------------------
327      !!                ***  ROUTINE lim_sbc_tau ***
328      !! 
329      !! ** Purpose : Update the ocean surface stresses due to the ice
330      !!         
331      !! ** Action  : * at each ice time step (every nn_fsbc time step):
332      !!                - compute the modulus of ice-ocean relative velocity
333      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
334      !!                      tmod_io = rhoco * | U_ice-U_oce |
335      !!                - update the modulus of stress at ocean surface
336      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
337      !!              * at each ocean time step (every kt):
338      !!                  compute linearized ice-ocean stresses as
339      !!                      Utau = tmod_io * | U_ice - pU_oce |
340      !!                using instantaneous current ocean velocity (usually before)
341      !!
342      !!    NB: - ice-ocean rotation angle no more allowed
343      !!        - here we make an approximation: taum is only computed every ice time step
344      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
345      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
346      !!
347      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
348      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
349      !!---------------------------------------------------------------------
350      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
351      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
352      !!
353      INTEGER  ::   ji, jj   ! dummy loop indices
354      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
355      REAL(wp) ::   zat_v, zvtau_ice, zv_t          !   -      -
356      !!---------------------------------------------------------------------
357      !
358      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
359!CDIR NOVERRCHK
360         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
361!CDIR NOVERRCHK
362            DO ji = fs_2, fs_jpim1
363               !                                               ! 2*(U_ice-U_oce) at T-point
364               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
365               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
366               !                                              ! |U_ice-U_oce|^2
367               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
368               !                                               ! update the ocean stress modulus
369               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt
370               tmod_io(ji,jj) = rhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
371            END DO
372         END DO
373         CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
374         !
375         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
376         vtau_oce(:,:) = vtau(:,:)
377         !
378      ENDIF
379      !
380      !                                      !==  every ocean time-step  ==!
381      !
382      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
383         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
384            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
385            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
386            !                                                   ! linearized quadratic drag formulation
387            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
388            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
389            !                                                   ! stresses at the ocean surface
390            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
391            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
392         END DO
393      END DO
394      CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
395      !
396      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
397         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
398     
399   END SUBROUTINE lim_sbc_tau
400
401
402   SUBROUTINE lim_sbc_init
403      !!-------------------------------------------------------------------
404      !!                  ***  ROUTINE lim_sbc_init  ***
405      !!             
406      !! ** Purpose : Preparation of the file ice_evolu for the output of
407      !!      the temporal evolution of key variables
408      !!
409      !! ** input   : Namelist namicedia
410      !!-------------------------------------------------------------------
411      REAL(wp) :: zsum, zarea
412      r1_rau0 = 1._wp / rau0
413      !
414      IF(lwp) WRITE(numout,*)
415      IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition'
416      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   '
417
418      !                                      ! allocate lim_sbc array
419      IF( lim_sbc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
420      !
421      r1_rdtice = 1. / rdt_ice
422      !
423      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
424      sice_0(:,:) = sice
425      !
426      IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea
427         WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
428            &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
429            soce_0(:,:) = 4._wp
430            sice_0(:,:) = 2._wp
431         END WHERE
432      ENDIF
433      ! clem modif
434      iatte(:,:) = 1._wp
435      oatte(:,:) = 1._wp
436      !
437      !                                      ! sea ice  with mass exchange
438         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
439         snwice_mass_b(:,:) = snwice_mass(:,:)
440         IF( .NOT. ln_rstart ) THEN           ! delete the initial ssh below sea-ice area
441         !
442         zarea     = glob_sum( e1e2t(:,:) )           ! interior global domain surface
443         zsum      = glob_sum( e1e2t(:,:) * ( snwice_mass(:,:) ) ) / zarea * r1_rau0
444         sshn(:,:) = sshn(:,:) - zsum 
445         sshb(:,:) = sshb(:,:) - zsum
446         ENDIF
447      !
448      !
449   END SUBROUTINE lim_sbc_init
450
451#else
452   !!----------------------------------------------------------------------
453   !!   Default option :        Dummy module       NO LIM 3.0 sea-ice model
454   !!----------------------------------------------------------------------
455CONTAINS
456   SUBROUTINE lim_sbc           ! Dummy routine
457   END SUBROUTINE lim_sbc
458#endif 
459
460   !!======================================================================
461END MODULE limsbc
Note: See TracBrowser for help on using the repository browser.