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

source: branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90 @ 1859

Last change on this file since 1859 was 1859, checked in by gm, 14 years ago

ticket:#665 step 2 & 3: heat content in qns & new forcing terms

  • Property svn:keywords set to Id
File size: 17.5 KB
Line 
1MODULE limsbc_2
2   !!======================================================================
3   !!                       ***  MODULE limsbc_2   ***
4   !!           computation of the flux at the sea ice/ocean interface
5   !!======================================================================
6   !! History :  1.0  !  2000-01  (H. Goosse) Original code
7   !!            2.0  !  2002-07  (C. Ethe, G. Madec) re-writing F90
8   !!             -   !  2006-07  (G. Madec) surface module
9   !!             -   !  2008-07  (C. Talandier,G.  Madec) 2D fields for soce and sice
10   !!            2.1  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
11   !!----------------------------------------------------------------------
12#if defined key_lim2
13   !!----------------------------------------------------------------------
14   !!   'key_lim2'                                    LIM 2.0 sea-ice model
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          ! ice   surface boundary condition
21   USE sbc_oce          ! ocean surface boundary condition
22   USE phycst           ! physical constants
23   USE albedo           ! albedo parameters
24   USE ice_2            ! LIM-2 sea-ice variables
25
26   USE lbclnk           ! ocean lateral boundary condition
27   USE in_out_manager   ! I/O manager
28   USE iom              !
29   USE prtctl           ! Print control
30   USE diaar5, ONLY :   lk_diaar5
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)  ::   epsi16 = 1.e-16  ! constant values
39   REAL(wp)  ::   rzero  = 0.e0   
40   REAL(wp)  ::   rone   = 1.e0
41   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r, sice_r   ! ocean and ice 2D constant salinity fields (used if lk_vvl=F)
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)
47   !! $Id$
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     : solar heat flux
65      !!              - qns     : non-solar heat flux including heat content of mass flux
66      !!              - emp     : mass flux
67      !!              - emps    : salt flux due to Freezing/Melting
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, INTENT(in) ::   kt    ! number of iteration
78      !!
79      INTEGER  ::   ji, jj                     ! dummy loop indices
80      INTEGER  ::   ifvt, idfr , iadv, i1mfr   ! local integers
81      INTEGER  ::   iflt, ifrdv, ial , ifral   !   -      -
82      INTEGER  ::   ii0, ii1, ij0, ij1         !   -      -
83      REAL(wp) ::   zqsr, zqns, zqhc, zemp     ! local scalars
84      REAL(wp) ::   zinda, zswitch, zcd        !   -      -
85      REAL(wp) ::   zfrldu, zutau, zu_io       !   -      -
86      REAL(wp) ::   zfrldv, zvtau, zv_io       !   -      -
87      REAL(wp) ::   zemp_snw, zfmm, zfsalt     !   -      -
88      REAL(wp) ::   zsang, zmod, zztmp, zfm
89      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp    ! 3D workspace
90      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! 2D workspace
91      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi, zqnsoce  !  -      -
92      !!---------------------------------------------------------------------
93           
94      IF( kt == nit000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
97         IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
98         !                              ! 2D fields for constant ice and ocean salinities
99         soce_r(:,:) = soce             !    in order to use different value in the Baltic sea
100         sice_r(:,:) = sice             !    which is much less salty than polar regions
101         !
102         IF( cp_cfg == "orca" ) THEN    ! ORCA configuration
103            IF( 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) ) = 0.50
107!!gm add here the R1 R05 and R025  cases
108!!          ELSEIF( jp_cfg == 1   ) THEN           !  ORCA_R1   configuration
109!!          ELSEIF( jp_cfg == 05  ) THEN           !  ORCA_R05  configuration
110!!          ELSEIF( jp_cfg == 025 ) THEN           !  ORCA_R025 configuration
111!!
112!!gm or better introduce the baltic change as a function of lat/lon of the baltic sea
113!!end gm
114            ENDIF
115         ENDIF
116         !
117      ENDIF
118
119      zqnsoce(:,:) = qns(:,:)      ! save non-solar flux prior to its modification by ice-ocean fluxes (diag.)
120      !
121      zswitch = 1       ! standard levitating sea-ice : salt exchange only
122      !
123!!gm ice embedment
124!      SELECT CASE( nn_ice_embd )       ! levitating/embedded sea-ice
125!      CASE( 0    )   ;   zswitch = 1       ! standard levitating sea-ice : salt exchange only
126!      CASE( 1, 2 )   ;   zswitch = 0       ! other levitating sea-ice or embedded sea-ice : salt and volume fluxes
127!      END SELECT                           !   
128!!gm end embedment
129
130      DO jj = 1, jpj
131         DO ji = 1, jpi
132            !                          !------------------------------------------!
133            !                          !      heat flux at the ocean surface      !
134            !                          !------------------------------------------!
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!!gm  attempt to understand and comment the tricky flags used here....
146!
147!gm      zinda   = 1.0 - AINT( pfrld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice adv, but before ice thermo)
148!gm      i1mfr   = 1.0 - AINT(  frld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice thermo)
149!
150!gm      IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda   ! = 1. if (snow and no ice at previous time) else 0. ???
151!gm      ELSE                             ;   ifvt = 0.      ! correspond to a overmelting of snow in surface ablation
152!gm      ENDIF                                               !
153!
154!gm      IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermo
155!gm      ELSE                                     ;   idfr = 1.   
156!gm      ENDIF
157!
158!!$      iflt    = zinda  * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous and pure ocean at current
159!
160!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
161!!$!                 snow no ice   ice         ice or nothing  lead fraction increases
162!!$!                 at previous   now           at previous
163!!$!                -> ice aera increases  ???         -> ice aera decreases ???
164!
165!!$            iadv    = ( 1  - i1mfr ) * zinda 
166!!$!                     pure ocean      ice at
167!!$!                     at current      previous
168!!$!                        -> = 1. if ice disapear between previous and current
169!
170!!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) ) 
171!!$!                            ice at     ???
172!!$!                            current         
173!!$!                         -> ???
174!
175!!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv
176!!$!                                                    ice disapear                           
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            !
186            ! - computation the non solar heat flux at ocean surface
187            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
188               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                                &
189               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdt_ice    &
190               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdt_ice
191
192            ! - store residual heat flux (put in the ocean at the next time-step)
193            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)   ! ???
194            !
195            ! - heat content of mass exchanged between ocean and sea-ice
196            zqhc = ( rdq_snw(ji,jj) + rdq_ice(ji,jj) ) * r1_rdt_ice    ! heat flux due to sown & ice heat content exchanges
197            !           
198            qsr(ji,jj) = zqsr                                          ! solar heat flux
199            qns(ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux
200 
201            !                          !------------------------------------------!
202            !                          !      mass flux at the ocean surface      !
203            !                          !------------------------------------------!
204            !
205            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
206#if defined key_coupled
207            !                                                       ! coupled mode:
208            zemp = + emp_tot(ji,jj)                              &       ! net mass flux over the grid cell (ice+ocean area)
209               &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )              ! minus the mass flux intercepted by sea-ice
210#else
211            !                                                       ! forced  mode:
212            zemp = + emp(ji,jj)     *         frld(ji,jj)      &         ! mass flux over open ocean fraction
213               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &         ! liquid precip. over ice reaches directly the ocean
214               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &         ! snow is intercepted by sea-ice (previous frld)
215#endif           
216            !
217            ! mass flux at the ocean/ice interface (sea ice fraction)
218            zemp_snw = rdm_snw(ji,jj) * r1_rdt_ice                  ! snow melting = pure water that enters the ocean
219            zfmm     = rdm_ice(ji,jj) * r1_rdt_ice                  ! Freezing minus Melting (F-M)
220
221            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
222            zfsalt = - sice_r(ji,jj) * zfmm                         ! F-M salt exchange
223            zcd    =   soce_r(ji,jj) * zfmm                         ! concentration/dilution term due to F-M
224            !
225            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux
226            ! salt and mass fluxes : non concentartion dilution term in salt flux  and add F-M term in volume flux
227            emps(ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange)
228            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (- F/M mass flux if no ice/ocean mass exchange)
229            !
230         END DO
231      END DO
232
233
234      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )     
235      CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )     
236      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) )
237
238      IF( lk_diaar5 ) THEN
239         CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * zrdtir )
240         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdm_ice(:,:) * zrdtir )
241         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdm_ice(:,:) * zrdtir )
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               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
276               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
277               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
278               zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) )
279               zfrldv = 0.5 * ( frld(ji,jj) + frld(ji  ,jj+1) )
280               ! ... update components of surface ocean stress (ice-cover wheighted)
281               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
282               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
283               ! ... module of ice-ocean stress at T-points (from I-point values)
284               zztmp = 0.25 * ( ztiomi(ji,jj) + ztiomi(ji+1,jj) + ztiomi(ji,jj+1) + ztiomi(ji+1,jj+1) )
285               ! ... update module of surface ocean stress (ice-cover wheighted)
286               taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1. - frld(ji,jj) ) * zztmp
287               !
288            END DO
289         END DO
290         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )         ! lateral boundary conditions
291         CALL lbc_lnk( taum, 'T',  1. )
292         !
293      ENDIF
294
295      !-----------------------------------------------!
296      !   Storing the transmitted variables           !
297      !-----------------------------------------------!
298!!gm where this is done ?????   ==>>> limthd_2   not logic ???
299!!gm      fr_i(:,:) = 1.0 - frld(:,:)       ! sea-ice fraction
300!!gm
301
302      IF ( lk_cpl ) THEN      ! coupled mode :
303         tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
304         !                                  ! snow/ice and ocean albedo
305         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb )
306         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
307         !
308         CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo
309      ENDIF
310
311      IF(ln_ctl) THEN
312         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
313         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
314         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
315            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
316         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
317      ENDIF 
318      !
319   END SUBROUTINE lim_sbc_2
320
321#else
322   !!----------------------------------------------------------------------
323   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
324   !!----------------------------------------------------------------------
325CONTAINS
326   SUBROUTINE lim_sbc_2         ! Dummy routine
327   END SUBROUTINE lim_sbc_2
328#endif 
329
330   !!======================================================================
331END MODULE limsbc_2
Note: See TracBrowser for help on using the repository browser.