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

source: branches/dev_1784_EVP/NEMO/LIM_SRC_2/limsbc_2.F90 @ 2300

Last change on this file since 2300 was 2119, checked in by cbricaud, 14 years ago

bugfix

  • Property svn:keywords set to Id
File size: 15.4 KB
Line 
1MODULE limsbc_2
2   !!======================================================================
3   !!                       ***  MODULE limsbc_2   ***
4   !!           computation of the flux at the sea ice/ocean interface
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   !!----------------------------------------------------------------------
11#if defined key_lim2
12   !!----------------------------------------------------------------------
13   !!   'key_lim2'                                    LIM 2.0 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_sbc_2  : flux at the ice / ocean interface
16   !!----------------------------------------------------------------------
17   USE par_oce          ! ocean parameters
18   USE dom_oce          ! ocean domain
19   USE sbc_ice          ! surface boundary condition: ice
20   USE sbc_oce          ! surface boundary condition: ocean
21   USE phycst           ! physical constants
22   USE ice_2            ! LIM2: ice variables
23
24   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
25   USE in_out_manager   ! I/O manager
26   USE diaar5, ONLY :   lk_diaar5
27   USE iom              ! IOM library
28   USE albedo           ! albedo parameters
29   USE prtctl           ! Print control
30   USE cpl_oasis3, ONLY : lk_cpl
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2
36
37   REAL(wp)  ::   r1_rdtice                    ! constant values
38   REAL(wp)  ::   epsi16 = 1.e-16              !     -      -
39   REAL(wp)  ::   rzero  = 0.e0                !     -      -
40   REAL(wp)  ::   rone   = 1.e0                !     -      -
41   !
42   REAL(wp), DIMENSION(jpi,jpj) ::   soce_r, sice_r   ! constant SSS and ice salinity used in levitating sea-ice case
43
44   !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE lim_sbc_2( kt )
54      !!-------------------------------------------------------------------
55      !!                ***  ROUTINE lim_sbc_2 ***
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      !!              - fr_i    : ice fraction
71      !!              - tn_ice  : sea-ice surface temperature
72      !!              - alb_ice : sea-ice alberdo (lk_cpl=T)
73      !!
74      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
75      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
76      !!---------------------------------------------------------------------
77      INTEGER ::   kt    ! number of iteration
78      !!
79      INTEGER  ::   ji, jj           ! dummy loop indices
80      INTEGER  ::   ii0, ii1, ij0, ij1         ! local integers
81      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
82      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
83      REAL(wp) ::   zqsr, zqns, zsang, zmod, zfm   ! local scalars
84      REAL(wp) ::   zinda, zfons, zemp, zztmp      !   -      -
85      REAL(wp) ::   zfrldu, zutau, zu_io           !   -      -
86      REAL(wp) ::   zfrldv, zvtau, zv_io           !   -      -
87      REAL(wp), DIMENSION(jpi,jpj)   ::   ztio_u, ztio_v    ! 2D workspace
88      REAL(wp), DIMENSION(jpi,jpj)   ::   ztiomi, zqnsoce   !  -     -
89      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp   ! 2D/3D workspace
90      !!---------------------------------------------------------------------
91     
92     
93      IF( kt == nit000 ) THEN
94         IF(lwp) WRITE(numout,*)
95         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
96         IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
97         !
98         r1_rdtice = 1. / rdt_ice
99         !
100         soce_r(:,:) = soce
101         sice_r(:,:) = sice
102         !
103         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN     !  ORCA_R2 configuration
104            ii0 = 145   ;   ii1 = 180        ! Baltic Sea
105            ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0
106                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 2.e0
107         ENDIF
108         !
109      ENDIF
110
111      !------------------------------------------!
112      !      heat flux at the ocean surface      !
113      !------------------------------------------!
114
115!!gm
116!!gm CAUTION   
117!!gm re-verifies the non solar expression, especially over open ocen
118!!gm
119      zqnsoce(:,:) = qns(:,:)
120      DO jj = 1, jpj
121         DO ji = 1, jpi
122            zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
123            ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
124            i1mfr   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) )    ) )
125            idfr    = 1.0   - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
126            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
127            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
128            iadv    = ( 1  - i1mfr ) * zinda
129            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
130            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
131
132!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time)
133!!$
134!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time)
135!!$
136!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ???
137!!$            ELSE                             ;   ifvt = 0.
138!!$            ENDIF
139!!$
140!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current
141!!$            ELSE                                     ;   idfr = 1.   
142!!$            ENDIF
143!!$
144!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current
145!!$
146!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
147!!$!                 snow no ice   ice         ice or nothing  lead fraction increases
148!!$!                 at previous   now           at previous
149!!$!                -> ice aera increases  ???         -> ice aera decreases ???
150!!$
151!!$            iadv    = ( 1  - i1mfr ) * zinda 
152!!$!                     pure ocean      ice at
153!!$!                     at current      previous
154!!$!                        -> = 1. if ice disapear between previous and current
155!!$
156!!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) ) 
157!!$!                            ice at     ???
158!!$!                            current         
159!!$!                         -> ???
160!!$
161!!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv
162!!$!                                                    ice disapear                           
163!!$
164!!$
165
166            ! solar flux at ocean surface
167#if defined key_coupled 
168            zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
169#else
170            zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
171#endif           
172            ! non solar heat flux at ocean surface
173            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
174               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                              &
175               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   &
176               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice
177            !
178            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
179            !
180            qsr  (ji,jj) = zqsr                                          ! solar heat flux
181            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
182         END DO
183      END DO
184
185      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:)                         )     
186      CALL iom_put( 'qns_io_cea'  , qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )     
187      CALL iom_put( 'qsr_io_cea'  , fstric(:,:) * ( 1.e0 - pfrld(:,:) )  )
188
189      !------------------------------------------!
190      !      mass flux at the ocean surface      !
191      !------------------------------------------!
192      DO jj = 1, jpj
193         DO ji = 1, jpi
194#if defined key_coupled
195            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode)
196            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )   &   ! atmosphere-ocean freshwater flux
197               &                  + rdmsnif(ji,jj) * r1_rdtice                   ! freshwater flux due to snow melting
198#else
199            ! freshwater exchanges at the ice-atmosphere / ocean interface (forced mode)
200            zemp = + emp(ji,jj)     *         frld(ji,jj)     &   ! e-p budget over open ocean fraction
201               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )   &   ! liquid precipitation reaches directly the ocean
202               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   ! (account for change in ice cover within the timestep
203               &   + rdmsnif(ji,jj) * r1_rdtice                   ! freshwaterflux due to snow melting
204#endif           
205            ! salt exchanges at the ice/ocean interface
206            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
207            !
208            ! convert the salt flux from ice into a freshwater flux from ocean
209            zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
210            !
211            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
212            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
213         END DO
214      END DO
215      !
216      IF( lk_diaar5 ) THEN
217         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice )
218         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * r1_rdtice )
219         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * r1_rdtice )
220      ENDIF
221
222      !------------------------------------------!
223      !    momentum flux at the ocean surface    !
224      !------------------------------------------!
225
226      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
227         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
228
229         ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
230!CDIR NOVERRCHK
231         DO jj = 1, jpj
232!CDIR NOVERRCHK
233            DO ji = 1, jpi
234               ! ... change the cosinus angle sign in the south hemisphere
235               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
236               ! ... ice velocity relative to the ocean at I-point
237               zu_io  = u_ice(ji,jj) - u_oce(ji,jj)
238               zv_io  = v_ice(ji,jj) - v_oce(ji,jj)
239               zmod   = SQRT( zu_io * zu_io + zv_io * zv_io )
240               zztmp  = rhoco * zmod
241               ! ... components of ice stress over ocean with a ice-ocean rotation angle (at I-point)
242               ztio_u(ji,jj) = zztmp * ( cangvg * zu_io - zsang * zv_io )
243               ztio_v(ji,jj) = zztmp * ( cangvg * zv_io + zsang * zu_io )
244               ! ... module of ice stress over ocean (at I-point)
245               ztiomi(ji,jj) = zztmp * zmod
246               !
247            END DO
248         END DO
249
250         DO jj = 2, jpjm1
251            DO ji = 2, jpim1   ! NO vector opt.
252               ! ... components of ice-ocean stress at U and V-points  (from I-point values)
253#if defined key_lim2_vp
254               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )      ! VP rheology
255               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
256#else
257               zutau  = ztio_u(ji,jj)                                      ! EVP rheology
258               zvtau  = ztio_v(ji,jj)
259#endif
260               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
261               zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) )
262               zfrldv = 0.5 * ( frld(ji,jj) + frld(ji  ,jj+1) )
263               ! ... update components of surface ocean stress (ice-cover wheighted)
264               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
265               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
266               ! ... module of ice-ocean stress at T-points (from I-point values)
267               zztmp = 0.25 * ( ztiomi(ji,jj) + ztiomi(ji+1,jj) + ztiomi(ji,jj+1) + ztiomi(ji+1,jj+1) )
268               ! ... update module of surface ocean stress (ice-cover wheighted)
269               taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1. - frld(ji,jj) ) * zztmp
270               !
271            END DO
272         END DO
273         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition
274         CALL lbc_lnk( taum, 'T',  1. )
275
276      ENDIF
277
278      IF( lk_cpl ) THEN           
279      !-----------------------------------------------!
280      !   Coupling variables                          !
281      !-----------------------------------------------!
282         tn_ice(:,:,1) = sist(:,:)           ! sea-ice surface temperature       
283         !                                   ! snow/ice and ocean albedo
284         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb )
285         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
286         CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo
287      ENDIF
288
289      IF(ln_ctl) THEN
290         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
291         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
292         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
293            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
294         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
295      ENDIF 
296      !
297   END SUBROUTINE lim_sbc_2
298
299#else
300   !!----------------------------------------------------------------------
301   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
302   !!----------------------------------------------------------------------
303CONTAINS
304   SUBROUTINE lim_sbc_2         ! Dummy routine
305   END SUBROUTINE lim_sbc_2
306#endif 
307
308   !!======================================================================
309END MODULE limsbc_2
Note: See TracBrowser for help on using the repository browser.