source: trunk/NEMO/LIM_SRC_3/limsbc.F90 @ 888

Last change on this file since 888 was 888, checked in by ctlod, 13 years ago

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113

File size: 19.3 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) ::   zqsr  , zqns     ! solar & non solar heat flux
80      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
81      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
82      REAL(wp) ::   zpme             ! freshwater exchanges at the ice/ocean interface
83      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
84      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
85      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
86      REAL(wp) ::   ztairx, ztairy   
87      REAL(wp) ::   zsfrldmx2, zsfrldmy2
88      REAL(wp) ::   zu_ice, zv_ice
89      REAL(wp) ::   ztglx , ztgly
90      REAL(wp) ::   ztaux , ztauy
91     
92#if defined key_coupled   
93      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky
94      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky
95#endif
96      REAL(wp) ::   zsang, zmod, zfm
97      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
98      REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes
99
100      !!---------------------------------------------------------------------
101     
102      IF( kt == nit000 ) THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'lim_sbc : LIM 3.0 sea-ice - surface boundary condition'
105         IF(lwp) WRITE(numout,*) '~~~~~~~   '
106      ENDIF
107
108      !------------------------------------------!
109      !      heat flux at the ocean surface      !
110      !------------------------------------------!
111      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
112      ! changed to old_frld and old ht_i
113       
114      DO jj = 1, jpj
115         DO ji = 1, jpi
116            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
117            ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif  (ji,jj) ) )  !subscripts are bad here
118            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( at_i(ji,jj)       ) ) )
119            idfr    = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) )
120            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
121            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
122            iadv    = ( 1  - i1mfr ) * zinda
123            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
124            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
125
126            ! switch --- 1.0 ---------------- 0.0 --------------------
127            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
128            ! zinda   | if pfrld = 1       | if pfrld < 1            |
129            !  -> ifvt| if pfrld old_ht_i
130            ! i1mfr   | if frld = 1        | if frld  < 1            |
131            ! idfr    | if frld <= pfrld    | if frld > pfrld        |
132            ! iflt    |
133            ! ial     |
134            ! iadv    |
135            ! ifral
136            ! ifrdv
137
138            !   computation the solar flux at ocean surface
139            zfcm1(ji,jj)   = pfrld(ji,jj) * qsr(ji,jj)  + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
140                ! fstric     Solar flux transmitted trough the ice
141                ! qsr        Net short wave heat flux on free ocean
142! new line
143            fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)
144
145            !  computation the non solar heat flux at ocean surface
146            zfcm2(ji,jj) = - zfcm1(ji,jj)                  &
147               &           + iflt    * ( fscmbq(ji,jj) )   & ! total abl -> fscmbq is given to the ocean
148! fscmbq and ffltbif are obsolete
149!              &           + iflt * ffltbif(ji,jj) !!! only if one category is used
150               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   &
151               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice                     &
152               &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!!
153               &           + fheat_rpo(ji,jj) & ! contribution from ridge formation
154               &           + fheat_res(ji,jj)
155                ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean
156                !         computed in limthd_zdf.F90
157                ! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t
158                ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok)
159                ! qldif   heat balance of the lead (or of the open ocean)
160                ! qfvbq   i think this is wrong!
161                ! ---> Array used to store energy in case of total lateral ablation
162                ! qfvbq latent heat uptake/release after accretion/ablation
163                ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
164
165            IF ( num_sal .EQ. 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + &
166                                  fhbri(ji,jj) ! new contribution due to brine drainage
167
168            ! bottom radiative component is sent to the computation of the
169            ! oceanic heat flux
170            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     
171
172            ! used to compute the oceanic heat flux at the next time step
173            qsr(ji,jj) = zfcm1(ji,jj)                                       ! solar heat flux
174            qns(ji,jj) = zfcm2(ji,jj) - fdtcn(ji,jj)                        ! non solar heat flux
175!                           ! fdtcn : turbulent oceanic heat flux
176
177
178            IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN
179               WRITE(numout,*) ' lim_sbc : heat fluxes '
180               WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx)
181               WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
182               WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx)
183               WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx)
184               WRITE(numout,*)
185               WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx)
186               WRITE(numout,*) ' zfcm2     : ', zfcm2(jiindx,jjindx)
187               WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx)
188               WRITE(numout,*) ' ifral     : ', ifral
189               WRITE(numout,*) ' ial       : ', ial 
190               WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx)
191               WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx)
192               WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) / rdt_ice
193               WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) / rdt_ice
194               WRITE(numout,*) ' ifrdv     : ', ifrdv
195               WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx)
196               WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx)
197               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) / rdt_ice
198               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) / rdt_ice
199               WRITE(numout,*) ' '
200               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
201               WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx)
202               WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx)
203               WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx)
204               WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
205            ENDIF
206         END DO
207      END DO
208       
209      !------------------------------------------!
210      !      mass flux at the ocean surface      !
211      !------------------------------------------!
212
213      DO jj = 1, jpj
214         DO ji = 1, jpi
215            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
216            !  -------------------------------------------------------------------------------------
217            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
218            !  Thus  FW  flux  =  External ( E-P+snow melt)
219            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
220            !                     Associated to Ice formation AND Ice melting
221            !                     Even if i see Ice melting as a FW and SALT flux
222            !       
223
224            !  computing freshwater exchanges at the ice/ocean interface
225            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj) )  &   !  evaporation over oceanic fraction
226               &   + tprecip(ji,jj) *         at_i(ji,jj)    &   !  total precipitation
227! old fashioned way               
228!              &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice
229               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   !  remov. snow precip over ice
230               &   - rdmsnif(ji,jj) / rdt_ice                &   !  freshwaterflux due to snow melting
231! new contribution from snow falling when ridging
232               &   + fmmec(ji,jj)
233           
234            !  computing salt exchanges at the ice/ocean interface
235            !  sice should be the same as computed with the ice model
236            zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
237! SOCE
238            zfons =  ( sss_m(ji,jj) - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
239           
240!CT useless            !  salt flux for constant salinity
241!CT useless            fsalt(ji,jj)      =  zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj)
242            !  salt flux for variable salinity
243            zinda             = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
244            !  correcting brine and salt fluxes
245            fsbri(ji,jj)      =  zinda*fsbri(ji,jj)
246            !  converting the salt fluxes from ice to a freshwater flux from ocean
247            fsalt_res(ji,jj)  =  fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 )
248            fseqv(ji,jj)      =  fseqv(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
249            fsbri(ji,jj)      =  fsbri(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
250            fsalt_rpo(ji,jj)  =  fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 )
251
252            !  freshwater mass exchange (positive to the ice, negative for the ocean ?)
253            !  actually it's a salt flux (so it's minus freshwater flux)
254            !  if sea ice grows, zfons is positive, fsalt also
255            !  POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN
256            !  POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1]
257
258            emp(ji,jj) = - zpme 
259         END DO
260      END DO
261
262      emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
263      IF (num_sal.eq.2) THEN
264         !In case of variable salinity the salt flux has to be accounted for differently
265         ! Brine drainage has to be added
266         emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
267      ENDIF
268     
269      IF( lk_dynspg_rl )    emp (:,:) = emps(:,:)      ! rigid-lid formulation : emp = emps
270
271      !------------------------------------------!
272      !    momentum flux at the ocean surface    !
273      !------------------------------------------!
274
275      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
276         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
277
278         ! ... ice stress over ocean with a ice-ocean rotation angle
279         DO jj = 2, jpjm1
280            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
281
282            DO ji = fs_2, fs_jpim1
283               ! computation of wind stress over ocean in X and Y direction
284               ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * utau(ji,jj)
285               ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * vtau(ji,jj)
286
287               zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
288               zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
289
290               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
291               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
292               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
293
294               ! quadratic drag formulation
295               ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
296               ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
297!
298!              ! IMPORTANT
299!              ! these lines are bound to prevent numerical oscillations
300!              ! in the ice-ocean stress
301!              ! They are physically ill-based. There is a cleaner solution
302!              ! to try (remember discussion in Paris Gurvan)
303!
304               ztglx   = ztglx * exp( - zmod / 0.5 )
305               ztgly   = ztglx * exp( - zmod / 0.5 )
306
307               tio_u(ji,jj) = ( ztairx + 1.0 * ztglx ) / 2.
308               tio_v(ji,jj) = ( ztairy + 1.0 * ztgly ) / 2.
309            END DO
310         END DO
311
312         DO jj = 1, jpjm1
313            DO ji = 1, fs_jpim1   ! vertor opt.
314               zfrldu = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points
315               zfrldv = MAX( freezn(ji,jj), freezn(ji+1,jj) )
316               ztaux = tio_u(ji,jj)
317               ztauy = tio_v(ji,jj)
318               utau(ji,jj) = (1.-zfrldu) * utau(ji,jj) + zfrldu * ztaux    ! stress at the ocean surface
319               vtau(ji,jj) = (1.-zfrldv) * vtau(ji,jj) + zfrldv * ztauy
320            END DO
321         END DO
322
323!new fashion         DO jj = 2, jpjm1
324!new fashion            DO ji = fs_2, fs_jpim1   ! vertor opt.
325!new fashion               ! ... ice-cover wheighted ice-ocean stress at U and V-points  (from I-point values)
326!new fashion               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
327!new fashion               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
328!new fashion               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
329!new fashion               zfrldu = 0.5 * ( frld (ji,jj) + frld (ji+1,jj  ) )
330!new fashion               zfrldv = 0.5 * ( frld (ji,jj) + frld (ji  ,jj+1) )
331!new fashion               ! update surface ocean stress
332!new fashion               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
333!new fashion               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
334!new fashion               !
335!new fashion            END DO
336!new fashion         END DO
337
338         ! boundary condition on the stress (utau,vtau)
339         CALL lbc_lnk( utau, 'U', -1. )
340         CALL lbc_lnk( vtau, 'V', -1. )
341
342      ENDIF
343
344      !-----------------------------------------------!
345      !   Storing the transmitted variables           !
346      !-----------------------------------------------!
347
348      freeze(:,:)   = at_i(:,:)             ! Sea ice cover           
349      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
350
351#if defined key_coupled           
352      !------------------------------------------------!
353      !    Computation of snow/ice and ocean albedo    !
354      !------------------------------------------------!
355      zalb  (:,:,:) = 0.e0
356      zalbp (:,:,:) = 0.e0
357
358      CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )
359
360      alb_ice(:,:,:) =  0.5 * zalbp(:,:,:) + 0.5 * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys)
361#endif
362
363      IF(ln_ctl) THEN
364         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
365         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
366         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
367            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
368         CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ')
369         CALL prt_ctl(tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl)
370      ENDIF
371   
372    END SUBROUTINE lim_sbc
373
374#else
375   !!----------------------------------------------------------------------
376   !!   Default option :        Dummy module       NO LIM 3.0 sea-ice model
377   !!----------------------------------------------------------------------
378CONTAINS
379   SUBROUTINE lim_sbc           ! Dummy routine
380   END SUBROUTINE lim_sbc
381#endif 
382
383   !!======================================================================
384END MODULE limsbc
Note: See TracBrowser for help on using the repository browser.