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 @ 917

Last change on this file since 917 was 917, checked in by rblod, 16 years ago

Correct a bug in revision 913 (ice/ocean stress computation)

File size: 17.7 KB
Line 
1MODULE limsbc
2   !!======================================================================
3   !!                       ***  MODULE limsbc   ***
4   !!           computation of the flux at the sea ice/ocean interface
5   !!======================================================================
6   !! History : 00-01 (H. Goosse) Original code
7   !!           02-07 (C. Ethe, G. Madec) re-writing F90
8   !!           06-07 (G. Madec) surface module
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                    LIM 3.0 sea-ice model
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   lim_sbc  : flux at the ice / ocean interface
16   !!----------------------------------------------------------------------
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 ocfzpt           ! surface ocean freezing point
24   USE ice_oce          ! sea-ice variable
25   USE ice              ! LIM sea-ice variables
26   USE iceini           ! ???
27   USE dynspg_oce       ! choice of the surface pressure gradient scheme
28
29   USE lbclnk           ! ocean lateral boundary condition
30   USE in_out_manager   ! I/O manager
31   USE albedo           ! albedo parameters
32   USE prtctl           ! Print control
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC lim_sbc       ! called by sbc_ice_lim
38
39   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
40   REAL(wp)  ::   rzero  = 0.e0   
41   REAL(wp)  ::   rone   = 1.e0
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2006)
47   !! $ Id: $
48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE lim_sbc( kt )
54      !!-------------------------------------------------------------------
55      !!                ***  ROUTINE lim_sbc ***
56      !! 
57      !! ** Purpose : Update surface ocean boundary condition over areas
58      !!      that are at least partially covered by sea-ice
59      !!         
60      !! ** Action  : - comput. of the momentum, heat and freshwater/salt
61      !!      fluxes at the ice-ocean interface.
62      !!              - Update
63      !!     
64      !! ** Outputs : - qsr    : sea heat flux:     solar
65      !!              - qns    : sea heat flux: non solar
66      !!              - emp    : freshwater budget: volume flux
67      !!              - emps   : freshwater budget: concentration/dillution
68      !!              - 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 iteration
75      !!
76      INTEGER  ::   ji, jj           ! dummy loop indices
77      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
78      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
79      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
80      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
81      REAL(wp) ::   zpme             ! freshwater exchanges at the ice/ocean interface
82      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
83      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
84      REAL(wp) ::   ztglx , ztgly    ! temporary scalar
85     
86#if defined key_coupled   
87      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky
88      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky
89#endif
90      REAL(wp) ::   zsang, zmod
91      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
92      REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes
93
94      !!---------------------------------------------------------------------
95     
96      IF( kt == nit000 ) THEN
97         IF(lwp) WRITE(numout,*)
98         IF(lwp) WRITE(numout,*) 'lim_sbc : LIM 3.0 sea-ice - surface boundary condition'
99         IF(lwp) WRITE(numout,*) '~~~~~~~   '
100      ENDIF
101
102      !------------------------------------------!
103      !      heat flux at the ocean surface      !
104      !------------------------------------------!
105      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
106      ! changed to old_frld and old ht_i
107       
108      DO jj = 1, jpj
109         DO ji = 1, jpi
110            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
111            ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif  (ji,jj) ) )  !subscripts are bad here
112            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( at_i(ji,jj)       ) ) )
113            idfr    = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) )
114            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
115            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
116            iadv    = ( 1  - i1mfr ) * zinda
117            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
118            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
119
120            ! switch --- 1.0 ---------------- 0.0 --------------------
121            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
122            ! zinda   | if pfrld = 1       | if pfrld < 1            |
123            !  -> ifvt| if pfrld old_ht_i
124            ! i1mfr   | if frld = 1        | if frld  < 1            |
125            ! idfr    | if frld <= pfrld    | if frld > pfrld        |
126            ! iflt    |
127            ! ial     |
128            ! iadv    |
129            ! ifral
130            ! ifrdv
131
132            !   computation the solar flux at ocean surface
133            zfcm1(ji,jj)   = pfrld(ji,jj) * qsr(ji,jj)  + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
134                ! fstric     Solar flux transmitted trough the ice
135                ! qsr        Net short wave heat flux on free ocean
136! new line
137            fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)
138
139            !  computation the non solar heat flux at ocean surface
140            zfcm2(ji,jj) = - zfcm1(ji,jj)                  &
141               &           + iflt    * ( fscmbq(ji,jj) )   & ! total abl -> fscmbq is given to the ocean
142! fscmbq and ffltbif are obsolete
143!              &           + iflt * ffltbif(ji,jj) !!! only if one category is used
144               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   &
145               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice                     &
146               &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!!
147               &           + fheat_rpo(ji,jj) & ! contribution from ridge formation
148               &           + fheat_res(ji,jj)
149                ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean
150                !         computed in limthd_zdf.F90
151                ! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t
152                ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok)
153                ! qldif   heat balance of the lead (or of the open ocean)
154                ! qfvbq   i think this is wrong!
155                ! ---> Array used to store energy in case of total lateral ablation
156                ! qfvbq latent heat uptake/release after accretion/ablation
157                ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
158
159            IF ( num_sal .EQ. 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + &
160                                  fhbri(ji,jj) ! new contribution due to brine drainage
161
162            ! bottom radiative component is sent to the computation of the
163            ! oceanic heat flux
164            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     
165
166            ! used to compute the oceanic heat flux at the next time step
167            qsr(ji,jj) = zfcm1(ji,jj)                                       ! solar heat flux
168            qns(ji,jj) = zfcm2(ji,jj) - fdtcn(ji,jj)                        ! non solar heat flux
169!                           ! fdtcn : turbulent oceanic heat flux
170
171
172            IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN
173               WRITE(numout,*) ' lim_sbc : heat fluxes '
174               WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx)
175               WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
176               WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx)
177               WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx)
178               WRITE(numout,*)
179               WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx)
180               WRITE(numout,*) ' zfcm2     : ', zfcm2(jiindx,jjindx)
181               WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
182               WRITE(numout,*) ' ifral     : ', ifral
183               WRITE(numout,*) ' ial       : ', ial 
184               WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx)
185               WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx)
186               WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) / rdt_ice
187               WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) / rdt_ice
188               WRITE(numout,*) ' ifrdv     : ', ifrdv
189               WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx)
190               WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx)
191               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) / rdt_ice
192               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) / rdt_ice
193               WRITE(numout,*) ' '
194               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
195               WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx)
196               WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx)
197               WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx)
198               WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
199            ENDIF
200         END DO
201      END DO
202       
203      !------------------------------------------!
204      !      mass flux at the ocean surface      !
205      !------------------------------------------!
206
207      DO jj = 1, jpj
208         DO ji = 1, jpi
209            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
210            !  -------------------------------------------------------------------------------------
211            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
212            !  Thus  FW  flux  =  External ( E-P+snow melt)
213            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
214            !                     Associated to Ice formation AND Ice melting
215            !                     Even if i see Ice melting as a FW and SALT flux
216            !       
217
218            !  computing freshwater exchanges at the ice/ocean interface
219            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj) )  &   !  evaporation over oceanic fraction
220               &   + tprecip(ji,jj) *         at_i(ji,jj)    &   !  total precipitation
221! old fashioned way               
222!              &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice
223               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   !  remov. snow precip over ice
224               &   - rdmsnif(ji,jj) / rdt_ice                &   !  freshwaterflux due to snow melting
225! new contribution from snow falling when ridging
226               &   + fmmec(ji,jj)
227           
228            !  computing salt exchanges at the ice/ocean interface
229            !  sice should be the same as computed with the ice model
230            zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
231! SOCE
232            zfons =  ( sss_m(ji,jj) - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
233           
234!CT useless            !  salt flux for constant salinity
235!CT useless            fsalt(ji,jj)      =  zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj)
236            !  salt flux for variable salinity
237            zinda             = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
238            !  correcting brine and salt fluxes
239            fsbri(ji,jj)      =  zinda*fsbri(ji,jj)
240            !  converting the salt fluxes from ice to a freshwater flux from ocean
241            fsalt_res(ji,jj)  =  fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 )
242            fseqv(ji,jj)      =  fseqv(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
243            fsbri(ji,jj)      =  fsbri(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
244            fsalt_rpo(ji,jj)  =  fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 )
245
246            !  freshwater mass exchange (positive to the ice, negative for the ocean ?)
247            !  actually it's a salt flux (so it's minus freshwater flux)
248            !  if sea ice grows, zfons is positive, fsalt also
249            !  POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN
250            !  POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1]
251
252            emp(ji,jj) = - zpme 
253         END DO
254      END DO
255
256      emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
257      IF (num_sal.eq.2) THEN
258         !In case of variable salinity the salt flux has to be accounted for differently
259         ! Brine drainage has to be added
260         emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
261      ENDIF
262     
263      IF( lk_dynspg_rl )    emp (:,:) = emps(:,:)      ! rigid-lid formulation : emp = emps
264
265      !------------------------------------------!
266      !    momentum flux at the ocean surface    !
267      !------------------------------------------!
268
269      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
270         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
271
272         ! ... ice stress over ocean with a ice-ocean rotation angle
273         DO jj = 1, jpj
274            ! ... change the sinus angle sign in the south hemisphere
275            zsang  = SIGN(1.e0, gphif(1,jj) ) * sangvg
276            DO ji = 1, jpi
277               ! ... ice velocity relative to the ocean
278               zu_io = u_ice(ji,jj) - u_oce(ji,jj)
279               zv_io = v_ice(ji,jj) - v_oce(ji,jj)
280               zmod  = SQRT( zu_io * zu_io + zv_io * zv_io ) 
281               ! quadratic drag formulation
282               ztglx = rhoco * zmod * ( cangvg * zu_io - zsang * zv_io )
283               ztgly = rhoco * zmod * ( cangvg * zv_io + zsang * zu_io )
284               ! IMPORTANT
285               ! these lines are bound to prevent numerical oscillations
286               ! in the ice-ocean stress
287               ! They are physically ill-based. There is a cleaner solution
288               ! to try (remember discussion in Paris Gurvan)
289               ztio_u(ji,jj) = ztglx * exp( - zmod / 0.5 ) 
290               ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 ) 
291               !
292            END DO
293         END DO
294
295         DO jj = 2, jpjm1
296            DO ji = fs_2, fs_jpim1   ! vertor opt.
297               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
298               zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) )
299               zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji  ,jj+1) )
300               ! update surface ocean stress
301               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj)
302               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj)
303               !
304            END DO
305         END DO
306
307         ! boundary condition on the stress (utau,vtau)
308         CALL lbc_lnk( utau, 'U', -1. )
309         CALL lbc_lnk( vtau, 'V', -1. )
310
311      ENDIF
312
313      !-----------------------------------------------!
314      !   Storing the transmitted variables           !
315      !-----------------------------------------------!
316
317      freeze(:,:)   = at_i(:,:)             ! Sea ice cover           
318      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
319
320#if defined key_coupled           
321      !------------------------------------------------!
322      !    Computation of snow/ice and ocean albedo    !
323      !------------------------------------------------!
324      zalb  (:,:,:) = 0.e0
325      zalbp (:,:,:) = 0.e0
326
327      CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )
328
329      alb_ice(:,:,:) =  0.5 * zalbp(:,:,:) + 0.5 * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys)
330#endif
331
332      IF(ln_ctl) THEN
333         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
334         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
335         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
336            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
337         CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ')
338         CALL prt_ctl(tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl)
339      ENDIF
340   
341    END SUBROUTINE lim_sbc
342
343#else
344   !!----------------------------------------------------------------------
345   !!   Default option :        Dummy module       NO LIM 3.0 sea-ice model
346   !!----------------------------------------------------------------------
347CONTAINS
348   SUBROUTINE lim_sbc           ! Dummy routine
349   END SUBROUTINE lim_sbc
350#endif 
351
352   !!======================================================================
353END MODULE limsbc
Note: See TracBrowser for help on using the repository browser.