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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/LIM_SRC_2/limsbc_2.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm

File size: 16.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 : 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_lim2
11   !!----------------------------------------------------------------------
12   !!   'key_lim2'                                    LIM 2.0 sea-ice model
13   !!----------------------------------------------------------------------
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
20   USE sbc_oce          ! surface boundary condition
21   USE phycst           ! physical constants
22   USE ice_2            ! LIM sea-ice variables
23
24   USE lbclnk           ! ocean lateral boundary condition
25   USE in_out_manager   ! I/O manager
26   USE diaar5, ONLY :   lk_diaar5
27   USE iom              !
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)  ::   epsi16 = 1.e-16  ! constant values
38   REAL(wp)  ::   rzero  = 0.e0   
39   REAL(wp)  ::   rone   = 1.e0
40   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r
41   REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)
47   !! $Id: limsbc_2.F90 1756 2009-11-25 14:15:20Z smasson $
48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
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  ::   ifvt, i1mfr, idfr               ! some switches
81      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
82      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
83      REAL(wp) ::   zrdtir           ! 1. / rdt_ice
84      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
85      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
86      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
87      REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface
88      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
89      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
90      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
91! interface 2D --> 3D
92      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky
93      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky
94      REAL(wp) ::   zsang, zmod, zztmp, zfm
95      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point
96      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point
97      REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model
98
99      !!---------------------------------------------------------------------
100     
101      zrdtir = 1. / rdt_ice
102     
103      IF( kt == nit000 ) THEN
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
106         IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
107
108         soce_r(:,:) = soce
109         sice_r(:,:) = sice
110         !
111         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
112            !                                        ! =======================
113            !                                        !  ORCA_R2 configuration
114            !                                        ! =======================
115            ii0 = 145   ;   ii1 = 180        ! Baltic Sea
116            ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0
117                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 2.e0
118         ENDIF
119         !
120      ENDIF
121
122      !------------------------------------------!
123      !      heat flux at the ocean surface      !
124      !------------------------------------------!
125
126!!gm
127!!gm CAUTION   
128!!gm re-verifies the non solar expression, especially over open ocen
129!!gm
130      zqnsoce(:,:) = qns(:,:)
131      DO jj = 1, jpj
132         DO ji = 1, jpi
133            zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
134            ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
135            i1mfr   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) )    ) )
136            idfr    = 1.0   - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
137            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
138            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
139            iadv    = ( 1  - i1mfr ) * zinda
140            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
141            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
142
143!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time)
144!!$
145!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time)
146!!$
147!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ???
148!!$            ELSE                             ;   ifvt = 0.
149!!$            ENDIF
150!!$
151!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current
152!!$            ELSE                                     ;   idfr = 1.   
153!!$            ENDIF
154!!$
155!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current
156!!$
157!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
158!!$!                 snow no ice   ice         ice or nothing  lead fraction increases
159!!$!                 at previous   now           at previous
160!!$!                -> ice aera increases  ???         -> ice aera decreases ???
161!!$
162!!$            iadv    = ( 1  - i1mfr ) * zinda 
163!!$!                     pure ocean      ice at
164!!$!                     at current      previous
165!!$!                        -> = 1. if ice disapear between previous and current
166!!$
167!!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) ) 
168!!$!                            ice at     ???
169!!$!                            current         
170!!$!                         -> ???
171!!$
172!!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv
173!!$!                                                    ice disapear                           
174!!$
175!!$
176
177            !   computation the solar flux at ocean surface
178#if defined key_coupled 
179            zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
180#else
181            zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
182#endif           
183            !  computation the non solar heat flux at ocean surface
184            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
185               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
186               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    &
187               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir
188
189            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
190           
191            qsr  (ji,jj) = zqsr                                          ! solar heat flux
192            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
193         END DO
194      END DO
195
196      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )     
197      CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )     
198      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) )
199
200      !------------------------------------------!
201      !      mass flux at the ocean surface      !
202      !------------------------------------------!
203
204!!gm
205!!gm CAUTION   
206!!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm
207!!gm
208      DO jj = 1, jpj
209         DO ji = 1, jpi
210           
211#if defined key_coupled
212          zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !
213             &   + rdmsnif(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting
214#else
215!!$            !  computing freshwater exchanges at the ice/ocean interface
216!!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction
217!!$               &   + tprecip(ji,jj)                           &   !  total precipitation
218!!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice
219!!$               &   - rdmsnif(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting
220            !  computing freshwater exchanges at the ice/ocean interface
221            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction
222               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
223               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
224               &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting
225               !                                                   !  ice-covered fraction:
226#endif           
227
228            !  computing salt exchanges at the ice/ocean interface
229            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir ) 
230           
231            !  converting the salt flux from ice to a freshwater flux from ocean
232            zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
233           
234            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
235            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
236
237         END DO
238      END DO
239
240      IF( lk_diaar5 ) THEN
241         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir )
242         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * zrdtir )
243         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * zrdtir )
244      ENDIF
245
246      !------------------------------------------!
247      !    momentum flux at the ocean surface    !
248      !------------------------------------------!
249
250      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
251         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
252
253         ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
254!CDIR NOVERRCHK
255         DO jj = 1, jpj
256!CDIR NOVERRCHK
257            DO ji = 1, jpi
258               ! ... change the cosinus angle sign in the south hemisphere
259               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
260               ! ... ice velocity relative to the ocean at I-point
261               zu_io  = u_ice(ji,jj) - u_oce(ji,jj)
262               zv_io  = v_ice(ji,jj) - v_oce(ji,jj)
263               zmod   = SQRT( zu_io * zu_io + zv_io * zv_io )
264               zztmp  = rhoco * zmod
265               ! ... components of ice stress over ocean with a ice-ocean rotation angle (at I-point)
266               ztio_u(ji,jj) = zztmp * ( cangvg * zu_io - zsang * zv_io )
267               ztio_v(ji,jj) = zztmp * ( cangvg * zv_io + zsang * zu_io )
268               ! ... module of ice stress over ocean (at I-point)
269               ztiomi(ji,jj) = zztmp * zmod
270               !
271            END DO
272         END DO
273
274         DO jj = 2, jpjm1
275            DO ji = 2, jpim1   ! NO vector opt.
276               ! ... components of ice-ocean stress at U and V-points  (from I-point values)
277               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
278               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
279               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
280               zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) )
281               zfrldv = 0.5 * ( frld(ji,jj) + frld(ji  ,jj+1) )
282               ! ... update components of surface ocean stress (ice-cover wheighted)
283               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
284               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
285               ! ... module of ice-ocean stress at T-points (from I-point values)
286               zztmp = 0.25 * ( ztiomi(ji,jj) + ztiomi(ji+1,jj) + ztiomi(ji,jj+1) + ztiomi(ji+1,jj+1) )
287               ! ... update module of surface ocean stress (ice-cover wheighted)
288               taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1. - frld(ji,jj) ) * zztmp
289               !
290            END DO
291         END DO
292
293         ! boundary condition on the stress (utau,vtau,taum)
294         CALL lbc_lnk( utau, 'U', -1. )
295         CALL lbc_lnk( vtau, 'V', -1. )
296         CALL lbc_lnk( taum, 'T',  1. )
297
298      ENDIF
299
300      !-----------------------------------------------!
301      !   Coupling variables                          !
302      !-----------------------------------------------!
303
304      IF ( lk_cpl ) THEN           
305         ! Ice surface temperature
306         tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
307         ! Computation of snow/ice and ocean albedo
308         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb )
309         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
310         CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo
311      ENDIF
312
313      IF(ln_ctl) THEN
314         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
315         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
316         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
317            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
318         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
319      ENDIF
320   
321    END SUBROUTINE lim_sbc_2
322
323#else
324   !!----------------------------------------------------------------------
325   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
326   !!----------------------------------------------------------------------
327CONTAINS
328   SUBROUTINE lim_sbc_2         ! Dummy routine
329   END SUBROUTINE lim_sbc_2
330#endif 
331
332   !!======================================================================
333END MODULE limsbc_2
Note: See TracBrowser for help on using the repository browser.