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

source: trunk/NEMO/LIM_SRC_3/limsbc.F90 @ 1684

Last change on this file since 1684 was 1684, checked in by rblod, 15 years ago

Correct a bug in LIM3 when averaging ocean surface velocity, see ticket #538

  • Property svn:keywords set to Id
File size: 26.4 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   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                    LIM 3.0 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_sbc  : flux at the ice / ocean interface
15   !!----------------------------------------------------------------------
16   USE oce
17   USE par_oce          ! ocean parameters
18   USE par_ice          ! ice parameters
19   USE dom_oce          ! ocean domain
20   USE sbc_ice          ! Surface boundary condition: sea-ice fields
21   USE sbc_oce          ! Surface boundary condition: ocean fields
22   USE phycst           ! physical constants
23   USE ice              ! LIM sea-ice variables
24   USE iceini           ! ???
25
26   USE lbclnk           ! ocean lateral boundary condition
27   USE in_out_manager   ! I/O manager
28   USE albedo           ! albedo parameters
29   USE prtctl           ! Print control
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   lim_sbc_flx   ! called by sbc_ice_lim
35   PUBLIC   lim_sbc_tau   ! called by sbc_ice_lim
36
37   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
38   REAL(wp)  ::   rzero  = 0.e0   
39   REAL(wp)  ::   rone   = 1.e0
40
41   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   !: air-ocean surface i- & j-stress              [N/m2]
42   REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              !: modulus of the ice-ocean relative velocity   [m/s]
43   REAL(wp), DIMENSION(jpi,jpj) ::   ssu_mb  , ssv_mb     !: before mean ocean surface currents           [m/s]
44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/LIM  3.2 ,  UCL-LOCEAN-IPSL (2009)
49   !! $Id$
50   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE lim_sbc_tau( kt, kcpl )
56      !!-------------------------------------------------------------------
57      !!                ***  ROUTINE lim_sbc_tau ***
58      !! 
59      !! ** Purpose : Update the ocean surface stresses due to the ice
60      !!         
61      !! ** Action  : - compute the ice-ocean stress depending on kcpl:
62      !!              fluxes at the ice-ocean interface.
63      !!     Case 0  :  old LIM-3 way, computed at ice time-step only
64      !!     Case 1  :  at each ocean time step the stresses are computed
65      !!                using the current ocean velocity (now)
66      !!     Case 2  :  combination of half case 0 + half case 1
67      !!     
68      !! ** Outputs : - utau   : sea surface i-stress (ocean referential)
69      !!              - vtau   : sea surface j-stress (ocean referential)
70      !!
71      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
72      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
73      !!---------------------------------------------------------------------
74      INTEGER ::   kt     ! number of ocean iteration
75      INTEGER ::   kcpl   ! ice-ocean coupling indicator: =0  LIM-3 old case
76      !                   !                               =1  stresses computed using now ocean velocity
77      !                   !                               =2  combination of 0 and 1 cases
78      !!
79      INTEGER  ::   ji, jj   ! dummy loop indices
80      REAL(wp) ::   zfrldu, zat_u, zu_ico, zutaui, zu_u, zu_ij, zu_im1j   ! temporary scalar
81      REAL(wp) ::   zfrldv, zat_v, zv_ico, zvtaui, zv_v, zv_ij, zv_ijm1   !    -         -
82      REAL(wp) ::   zsang                                                 !    -         -
83      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
84#if defined key_coupled   
85      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky
86      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky
87#endif
88     !!---------------------------------------------------------------------
89
90      IF( kt == nit000 ) THEN
91         IF(lwp) WRITE(numout,*)
92         IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM 3.0 sea-ice - surface ocean momentum fluxes'
93         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
94      ENDIF
95
96      SELECT CASE( kcpl )
97         !                                        !--------------------------------!
98      CASE( 0 )                                   !  LIM 3 old stress computation  !  (at ice timestep only)
99         !                                        !--------------------------------!
100         DO jj = 2, jpjm1                             !* modulus of the ice-ocean velocity
101            DO ji = fs_2, fs_jpim1
102               zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j)
103               zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j)
104               zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  )
105               zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1)
106               tmod_io(ji,jj) = SQRT(  0.5 * ( zu_ij * zu_ij + zu_im1j * zu_im1j   &
107                  &                          + zv_ij * zv_ij + zv_ijm1 * zv_ijm1  ) )
108            END DO
109         END DO
110         CALL lbc_lnk( tmod_io, 'T', 1. )   ! lateral boundary condition
111         !
112         !                                             !* ice stress over ocean with a ice-ocean rotation angle
113         DO jj = 1, jpjm1
114            zsang  = SIGN( 1.e0, gphif(1,jj) ) * sangvg         ! change the sinus angle sign in the south hemisphere
115            DO ji = 1, fs_jpim1
116               zu_u = u_ice(ji,jj) - u_oce(ji,jj)               ! ice velocity relative to the ocean
117               zv_v = v_ice(ji,jj) - v_oce(ji,jj)
118               !                                                ! quadratic drag formulation with rotation
119!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1)
120               zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v ) 
121               zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u ) 
122               !                                                ! bound for too large stress values
123               ! IMPORTANT: the exponential below prevents numerical oscillations in the ice-ocean stress
124               ! This is not physically based. A cleaner solution is offer in CASE kcpl=2
125               ztio_u(ji,jj) = zutaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji+1,jj) )  )
126               ztio_v(ji,jj) = zvtaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji,jj+1) )  ) 
127            END DO
128         END DO
129         DO jj = 2, jpjm1                                       ! stress at the surface of the ocean
130            DO ji = fs_2, fs_jpim1   ! vertor opt.
131               zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) )   ! open-ocean fraction at U- & V-points (from T-point values)
132               zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji  ,jj+1) )
133               !                                                    ! update surface ocean stress
134               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj)
135               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj)
136            END DO
137         END DO
138         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
139
140         !
141         !                                        !--------------------------------!
142      CASE( 1 )                                   !  Use the now velocity          !  (at each ocean timestep)
143         !                                        !--------------------------------!
144         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
145            utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step
146            vtau_oce(:,:) = vtau(:,:)
147            DO jj = 2, jpjm1                          !* modulus of the ice-ocean velocity
148               DO ji = fs_2, fs_jpim1
149                  zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j)
150                  zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j)
151                  zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  )
152                  zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1)
153                  tmod_io(ji,jj) = SQRT( 0.5 * (  zu_ij * zu_ij + zu_im1j * zu_im1j   &
154                     &                          + zv_ij * zv_ij + zv_ijm1 * zv_ijm1 ) )
155               END DO
156            END DO
157            CALL lbc_lnk( tmod_io, 'T', 1. )          ! lateral boundary condition
158         ENDIF
159         !
160        DO jj = 2, jpjm1                              !* ice stress over ocean with a ice-ocean rotation angle
161            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg        ! change the sinus angle sign in the south hemisphere
162            DO ji = fs_2, fs_jpim1
163               zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5   ! ice area at u and V-points
164               zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5
165               !                                                ! (u,v) ice-ocean velocity at (U,V)-point, resp.
166               zu_u   = u_ice(ji,jj) - ub(ji,jj,1)
167               zv_v   = v_ice(ji,jj) - vb(ji,jj,1)
168               !                                                ! quadratic drag formulation with rotation
169!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1)
170               zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v ) 
171               zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u ) 
172               !                                                   ! stress at the ocean surface
173               utau(ji,jj) = ( 1.- zat_u ) * utau_oce(ji,jj) + zat_u * zutaui
174               vtau(ji,jj) = ( 1.- zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui
175            END DO
176         END DO
177         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
178
179         !
180         !                                        !--------------------------------!
181      CASE( 2 )                                   !  mixed 0 and 2 cases           !  (at each ocean timestep)
182         !                                        !--------------------------------!
183         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
184            utau_oce(:,:) = utau (:,:)               !* save the air-ocean stresses at ice time-step
185            vtau_oce(:,:) = vtau (:,:)
186            ssu_mb  (:,:) = ssu_m(:,:)               !* save the ice-ocean velocity at ice time-step
187            ssv_mb  (:,:) = ssv_m(:,:)
188            DO jj = 2, jpjm1                         !* modulus of the ice-ocean velocity
189               DO ji = fs_2, fs_jpim1
190                  zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j)
191                  zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j)
192                  zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  )
193                  zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1)
194                  tmod_io(ji,jj) = SQRT( 0.5 * ( zu_ij * zu_ij + zu_im1j * zu_im1j   &
195                     &                         + zv_ij * zv_ij + zv_ijm1 * zv_ijm1 ) )
196               END DO
197            END DO
198            CALL lbc_lnk( tmod_io, 'T', 1. )
199         ENDIF
200         DO jj = 2, jpjm1                            !* ice stress over ocean with a ice-ocean rotation angle
201            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
202            DO ji = fs_2, fs_jpim1
203               zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5     ! ice area at u and V-points
204               zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 
205               !
206               zu_ico = u_ice(ji,jj) - 0.5 * ( ub(ji,jj,1) + ssu_mb(ji,jj) )   ! ice-oce velocity using ub and ssu_mb
207               zv_ico = v_ice(ji,jj) - 0.5 * ( vb(ji,jj,1) + ssv_mb(ji,jj) )
208               !                                        ! quadratic drag formulation with rotation
209!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1)
210               zutaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_ico - zsang * zv_ico )
211               zvtaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_ico + zsang * zu_ico )
212               !
213               utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface
214               vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui
215            END DO
216         END DO
217         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
218         !
219      END SELECT
220
221      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
222         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
223     
224   END SUBROUTINE lim_sbc_tau
225
226
227   SUBROUTINE lim_sbc_flx( kt )
228      !!-------------------------------------------------------------------
229      !!                ***  ROUTINE lim_sbc_flx ***
230      !! 
231      !! ** Purpose :   Update the surface ocean boundary condition for heat
232      !!              salt and mass over areas where sea-ice is non-zero
233      !!         
234      !! ** Action  : - computes the heat and freshwater/salt fluxes
235      !!              at the ice-ocean interface.
236      !!              - Update the ocean sbc
237      !!     
238      !! ** Outputs : - qsr     : sea heat flux:     solar
239      !!              - qns     : sea heat flux: non solar
240      !!              - emp     : freshwater budget: volume flux
241      !!              - emps    : freshwater budget: concentration/dillution
242      !!              - fr_i    : ice fraction
243      !!              - tn_ice  : sea-ice surface temperature
244      !!              - alb_ice : sea-ice alberdo (lk_cpl=T)
245      !!
246      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
247      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
248      !!---------------------------------------------------------------------
249      INTEGER ::   kt    ! number of iteration
250      !!
251      INTEGER  ::   ji, jj           ! dummy loop indices
252      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
253      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
254      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
255      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
256      REAL(wp) ::   zpme             ! freshwater exchanges at the ice/ocean interface
257      REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes
258#if defined key_coupled   
259      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky
260      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky
261#endif
262      !!---------------------------------------------------------------------
263
264      IF( kt == nit000 ) THEN
265         IF(lwp) WRITE(numout,*)
266         IF(lwp) WRITE(numout,*) 'lim_sbc_flx : LIM 3.0 sea-ice - heat salt and mass ocean surface fluxes'
267         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
268      ENDIF
269
270      !------------------------------------------!
271      !      heat flux at the ocean surface      !
272      !------------------------------------------!
273      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
274      ! changed to old_frld and old ht_i
275
276      DO jj = 1, jpj
277         DO ji = 1, jpi
278            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
279            ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif  (ji,jj) ) )  !subscripts are bad here
280            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( at_i(ji,jj)       ) ) )
281            idfr    = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) )
282            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
283            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
284            iadv    = ( 1  - i1mfr ) * zinda
285            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
286            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
287
288            ! switch --- 1.0 ---------------- 0.0 --------------------
289            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
290            ! zinda   | if pfrld = 1       | if pfrld < 1            |
291            !  -> ifvt| if pfrld old_ht_i
292            ! i1mfr   | if frld = 1        | if frld  < 1            |
293            ! idfr    | if frld <= pfrld    | if frld > pfrld        |
294            ! iflt    |
295            ! ial     |
296            ! iadv    |
297            ! ifral
298            ! ifrdv
299
300            !   computation the solar flux at ocean surface
301            zfcm1(ji,jj)   = pfrld(ji,jj) * qsr(ji,jj)  + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
302            ! fstric     Solar flux transmitted trough the ice
303            ! qsr        Net short wave heat flux on free ocean
304            ! new line
305            fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)
306
307            !  computation the non solar heat flux at ocean surface
308            zfcm2(ji,jj) = - zfcm1(ji,jj)                  &
309               &           + iflt    * ( fscmbq(ji,jj) )   & ! total abl -> fscmbq is given to the ocean
310               ! fscmbq and ffltbif are obsolete
311               !              &           + iflt * ffltbif(ji,jj) !!! only if one category is used
312               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   &
313               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice                     &
314               &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!!
315               &           + fheat_rpo(ji,jj) & ! contribution from ridge formation
316               &           + fheat_res(ji,jj)
317            ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean
318            !         computed in limthd_zdf.F90
319            ! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t
320            ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok)
321            ! qldif   heat balance of the lead (or of the open ocean)
322            ! qfvbq   i think this is wrong!
323            ! ---> Array used to store energy in case of total lateral ablation
324            ! qfvbq latent heat uptake/release after accretion/ablation
325            ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
326
327            IF ( num_sal .EQ. 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + &
328               fhbri(ji,jj) ! new contribution due to brine drainage
329
330            ! bottom radiative component is sent to the computation of the
331            ! oceanic heat flux
332            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     
333
334            ! used to compute the oceanic heat flux at the next time step
335            qsr(ji,jj) = zfcm1(ji,jj)                                       ! solar heat flux
336            qns(ji,jj) = zfcm2(ji,jj) - fdtcn(ji,jj)                        ! non solar heat flux
337            !                           ! fdtcn : turbulent oceanic heat flux
338
339            !!gm   this IF prevents the vertorisation of the whole loop
340            IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN
341               WRITE(numout,*) ' lim_sbc : heat fluxes '
342               WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx)
343               WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
344               WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx)
345               WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx)
346               WRITE(numout,*)
347               WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx)
348               WRITE(numout,*) ' zfcm2     : ', zfcm2(jiindx,jjindx)
349               WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
350               WRITE(numout,*) ' ifral     : ', ifral
351               WRITE(numout,*) ' ial       : ', ial 
352               WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx)
353               WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx)
354               WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) / rdt_ice
355               WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) / rdt_ice
356               WRITE(numout,*) ' ifrdv     : ', ifrdv
357               WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx)
358               WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx)
359               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) / rdt_ice
360               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) / rdt_ice
361               WRITE(numout,*) ' '
362               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
363               WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx)
364               WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx)
365               WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx)
366               WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
367            ENDIF
368            !!gm   end
369         END DO
370      END DO
371
372      !------------------------------------------!
373      !      mass flux at the ocean surface      !
374      !------------------------------------------!
375
376!!gm   optimisation: this loop have to be merged with the previous one
377      DO jj = 1, jpj
378         DO ji = 1, jpi
379            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
380            !  -------------------------------------------------------------------------------------
381            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
382            !  Thus  FW  flux  =  External ( E-P+snow melt)
383            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
384            !                     Associated to Ice formation AND Ice melting
385            !                     Even if i see Ice melting as a FW and SALT flux
386            !       
387
388            !  computing freshwater exchanges at the ice/ocean interface
389            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj) )  &   !  evaporation over oceanic fraction
390               &   + tprecip(ji,jj) *         at_i(ji,jj)    &   !  total precipitation
391               ! old fashioned way               
392               !              &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice
393               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   !  remov. snow precip over ice
394               &   - rdmsnif(ji,jj) / rdt_ice                &   !  freshwaterflux due to snow melting
395               ! new contribution from snow falling when ridging
396               &   + fmmec(ji,jj)
397
398            !  computing salt exchanges at the ice/ocean interface
399            !  sice should be the same as computed with the ice model
400            zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
401            ! SOCE
402            zfons =  ( sss_m(ji,jj) - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
403
404            !CT useless            !  salt flux for constant salinity
405            !CT useless            fsalt(ji,jj)      =  zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj)
406            !  salt flux for variable salinity
407            zinda             = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
408            !  correcting brine and salt fluxes
409            fsbri(ji,jj)      =  zinda*fsbri(ji,jj)
410            !  converting the salt fluxes from ice to a freshwater flux from ocean
411            fsalt_res(ji,jj)  =  fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 )
412            fseqv(ji,jj)      =  fseqv(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
413            fsbri(ji,jj)      =  fsbri(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
414            fsalt_rpo(ji,jj)  =  fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 )
415
416            !  freshwater mass exchange (positive to the ice, negative for the ocean ?)
417            !  actually it's a salt flux (so it's minus freshwater flux)
418            !  if sea ice grows, zfons is positive, fsalt also
419            !  POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN
420            !  POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1]
421
422            emp(ji,jj) = - zpme 
423         END DO
424      END DO
425
426      IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux
427         emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
428      ELSE                         ! constant ice salinity:
429         emps(:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
430      ENDIF
431
432      !-----------------------------------------------!
433      !   Storing the transmitted variables           !
434      !-----------------------------------------------!
435
436      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
437      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
438
439#if defined key_coupled           
440      !------------------------------------------------!
441      !    Computation of snow/ice and ocean albedo    !
442      !------------------------------------------------!
443      zalb  (:,:,:) = 0.e0
444      zalbp (:,:,:) = 0.e0
445
446      CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )
447
448      alb_ice(:,:,:) =  0.5 * zalbp(:,:,:) + 0.5 * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys)
449#endif
450
451      IF(ln_ctl) THEN
452         CALL prt_ctl( tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns , clinfo2=' qns     : ' )
453         CALL prt_ctl( tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps, clinfo2=' emps    : ' )
454         CALL prt_ctl( tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ' )
455         CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
456      ENDIF
457      !
458   END SUBROUTINE lim_sbc_flx
459
460#else
461   !!----------------------------------------------------------------------
462   !!   Default option :        Dummy module       NO LIM 3.0 sea-ice model
463   !!----------------------------------------------------------------------
464CONTAINS
465   SUBROUTINE lim_sbc           ! Dummy routine
466   END SUBROUTINE lim_sbc
467#endif 
468
469   !!======================================================================
470END MODULE limsbc
Note: See TracBrowser for help on using the repository browser.