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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90 @ 2319

Last change on this file since 2319 was 2319, checked in by cbricaud, 13 years ago

put new EVP rheology lost during the merge

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