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.F90 in trunk/NEMO/LIM_SRC – NEMO

source: trunk/NEMO/LIM_SRC/limsbc.F90 @ 702

Last change on this file since 702 was 702, checked in by smasson, 17 years ago

add first set of new surface module, see ticket:3

  • Property svn:executable set to *
File size: 11.9 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_ice_lim
11   !!----------------------------------------------------------------------
12   !!   'key_ice_lim'                                     LIM sea-ice model
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   lim_sbc  : 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 ocfzpt           ! surface ocean freezing point
23   USE ice_oce          ! sea-ice variable
24   USE ice              ! LIM sea-ice variables
25   USE iceini           ! ???
26   USE dynspg_oce       ! choice of the surface pressure gradient scheme
27
28   USE lbclnk           ! ocean lateral boundary condition
29   USE in_out_manager   ! I/O manager
30   USE albedo           ! albedo parameters
31   USE prtctl           ! Print control
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC lim_sbc       ! called by lim_step
37
38   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
39   REAL(wp)  ::   rzero  = 0.e0   
40   REAL(wp)  ::   rone   = 1.e0
41
42   !! * Substitutions
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)
46   !! $Header: $
47   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49
50CONTAINS
51
52   SUBROUTINE lim_sbc( kt )
53      !!-------------------------------------------------------------------
54      !!                ***  ROUTINE lim_sbc ***
55      !! 
56      !! ** Purpose : Update surface ocean boundary condition over areas
57      !!      that are at least partially covered by sea-ice
58      !!         
59      !! ** Action  : - comput. of the momentum, heat and freshwater/salt
60      !!      fluxes at the ice-ocean interface.
61      !!              - Update
62      !!     
63      !! ** Outputs : - fsolar  : solar heat flux at sea ice/ocean interface
64      !!              - fnsolar : non solar heat flux
65      !!              - fsalt   : salt flux at sea ice/ocean interface
66      !!              - fmass   : freshwater flux at sea ice/ocean interface
67      !!
68      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
69      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
70      !!---------------------------------------------------------------------
71      INTEGER ::   kt    ! number of iteration
72      !!
73      INTEGER  ::   ji, jj           ! dummy loop indices
74      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
75      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
76      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
77      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
78      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
79      REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface
80      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
81      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
82      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
83#if defined key_coupled   
84      REAL(wp), DIMENSION(jpi,jpj) ::   zalb     ! albedo of ice under overcast sky
85      REAL(wp), DIMENSION(jpi,jpj) ::   zalcn    ! albedo of ocean under overcast sky
86      REAL(wp), DIMENSION(jpi,jpj) ::   zalbp    ! albedo of ice under clear sky
87      REAL(wp), DIMENSION(jpi,jpj) ::   zaldum   ! albedo of ocean under clear sky
88#endif
89      REAL(wp) ::   zsang, zmod, zfm
90      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
91
92      !!---------------------------------------------------------------------
93     
94      IF( kt == nit000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'lim_sbc : LIM sea-ice - surface boundary condition'
97         IF(lwp) WRITE(numout,*) '~~~~~~~   '
98      ENDIF
99
100      !------------------------------------------!
101      !      heat flux at the ocean surface      !
102      !------------------------------------------!
103
104!!gm
105!!gm CAUTION   
106!!gm re-verifies the non solar expression, especially over open ocen
107!!gm
108      DO jj = 1, jpj
109         DO ji = 1, jpi
110            zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
111            ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
112            i1mfr   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) )    ) )
113            idfr    = 1.0   - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
114            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
115            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
116            iadv    = ( 1  - i1mfr ) * zinda
117            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
118            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
119            !   computation the solar flux at ocean surface
120            zqsr    = pfrld(ji,jj) * qsr(ji,jj)  + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
121            !  computation the non solar heat flux at ocean surface
122            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
123               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
124               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   &
125               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice
126
127            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
128           
129            qsr  (ji,jj) = zqsr                                          ! solar heat flux
130            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
131         END DO
132      END DO
133 
134      !------------------------------------------!
135      !      mass flux at the ocean surface      !
136      !------------------------------------------!
137
138!!gm
139!!gm CAUTION   
140!!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm
141!!gm
142      DO jj = 1, jpj
143         DO ji = 1, jpi
144           
145            !  computing freshwater exchanges at the ice/ocean interface
146            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction
147               !                                                   !  ice-covered fraction:
148               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
149               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
150               &   + rdmsnif(ji,jj) / rdt_ice                      !  freshwaterflux due to snow melting
151
152!!gm old    zemp = + evap(ji,jj) * frld(ji,jj)                 &   !  evaporation over oceanic fraction
153!!gm old       &   - tprecip(ji,jj)                            &   !  total precipitation
154!!gm old       &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  remov. snow precip over ice (store over sea-ice)
155!!gm old       &   + rdmsnif(ji,jj) / rdt_ice                      !  freshwaterflux due to snow melting
156           
157            !  computing salt exchanges at the ice/ocean interface
158            zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
159           
160            !  converting the salt flux from ice to a freshwater flux from ocean
161!!gm old    zfm  = zfons / ( sss_io(ji,jj) + epsi16 )
162            zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
163           
164            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
165            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
166
167         END DO
168      END DO
169
170      IF( lk_dynspg_rl )    emp (:,:) = emps(:,:)      ! rigid-lid formulation : emp = emps
171
172      !------------------------------------------!
173      !    momentum flux at the ocean surface    !
174      !------------------------------------------!
175
176      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
177         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
178
179         ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
180!CDIR NOVERRCHK
181         DO jj = 1, jpj
182!CDIR NOVERRCHK
183            DO ji = 1, jpi
184               ! ... change the cosinus angle sign in the south hemisphere
185               zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg 
186               ! ... ice velocity relative to the ocean
187               zu_io  = u_ice(ji,jj) - u_oce(ji,jj)
188               zv_io  = v_ice(ji,jj) - v_oce(ji,jj)
189               zmod   = SQRT( zu_io * zu_io + zv_io * zv_io )
190               ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
191               ztio_u(ji,jj) = cw * zmod * ( cangvg * zu_io - zsang * zv_io )
192               ztio_v(ji,jj) = cw * zmod * ( cangvg * zv_io + zsang * zu_io )
193               !
194            END DO
195         END DO
196
197         DO jj = 2, jpjm1
198            DO ji = 2, jpim1
199               ! ... ice-cover wheighted ice-ocean stress at U and V-points  (from I-point values)
200               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
201               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_u(ji+1,jj+1) )
202               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
203               zfrldu = 0.5 * ( frld (ji,jj) + frld (ji+1,jj  ) )
204               zfrldv = 0.5 * ( frld (ji,jj) + frld (ji  ,jj+1) )
205               ! update surface ocean stress
206               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
207               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
208               !
209            END DO
210         END DO
211
212         ! boundary condition on the stress (utau,vtau)
213         CALL lbc_lnk( utau, 'U', -1. )
214         CALL lbc_lnk( vtau, 'V', -1. )
215
216      ENDIF
217
218#if defined key_coupled           
219      !------------------------------------------------!
220      !    Computation of snow/ice and ocean albedo    !
221      !------------------------------------------------!
222      zalb  (:,:) = 0.e0
223      zalcn (:,:) = 0.e0
224      zalbp (:,:) = 0.e0
225      zaldum(:,:) = 0.e0
226
227      CALL flx_blk_albedo( zalb, zalcn, zalbp, zaldum )
228
229      alb_ice(:,:) =  0.5 * zalbp(:,:) + 0.5 * zalb (:,:)   ! Ice albedo (mean clear and overcast skys)
230#endif
231
232      IF(ln_ctl) THEN
233         CALL prt_ctl(tab2d_1=qsr , clinfo1=' lim_sbc: qsr  : ', tab2d_2=qns , clinfo2=' qns  : ')
234         CALL prt_ctl(tab2d_1=emp , clinfo1=' lim_sbc: emp  : ', tab2d_2=emps, clinfo2=' emps : ')
235         CALL prt_ctl(tab2d_1=utau, clinfo1=' lim_sbc: utau : ', tab2d_2=vtau, clinfo2=' vtau : ')
236!!gm tn_ice ???  what is it now?
237!!gm     CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ', tab2d_2=tn_ice , clinfo2=' tn_ice  : ')
238      ENDIF
239   
240    END SUBROUTINE lim_sbc
241
242#else
243   !!----------------------------------------------------------------------
244   !!   Default option :        Dummy module           NO LIM sea-ice model
245   !!----------------------------------------------------------------------
246CONTAINS
247   SUBROUTINE lim_sbc         ! Dummy routine
248   END SUBROUTINE lim_sbc
249#endif 
250
251   !!======================================================================
252END MODULE limsbc
Note: See TracBrowser for help on using the repository browser.