Changeset 918 for trunk/NEMO/LIM_SRC_3/limsbc.F90
- Timestamp:
- 2008-05-09T12:00:09+02:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limsbc.F90
r917 r918 4 4 !! computation of the flux at the sea ice/ocean interface 5 5 !!====================================================================== 6 !! History : 00-01 (H. Goosse) Original code7 !! 02-07 (C. Ethe, G. Madec) re-writing F908 !! 06-07 (G. Madec) surface module6 !! History : - ! 2006-07 (M. Vancoppelle) LIM3 original code 7 !! 3.0 ! 2008-03 (C. Tallandier) surface module 8 !! - ! 2008-04 (C. Tallandier) split in 2 + new ice-ocean coupling 9 9 !!---------------------------------------------------------------------- 10 10 #if defined key_lim3 … … 12 12 !! 'key_lim3' LIM 3.0 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 !!----------------------------------------------------------------------15 14 !! lim_sbc : flux at the ice / ocean interface 16 15 !!---------------------------------------------------------------------- 16 USE oce 17 17 USE par_oce ! ocean parameters 18 18 USE par_ice ! ice parameters … … 35 35 PRIVATE 36 36 37 PUBLIC lim_sbc ! called by sbc_ice_lim 37 PUBLIC lim_sbc_flx ! called by sbc_ice_lim 38 PUBLIC lim_sbc_tau ! called by sbc_ice_lim 38 39 39 40 REAL(wp) :: epsi16 = 1.e-16 ! constant values … … 41 42 REAL(wp) :: rone = 1.e0 42 43 44 REAL(wp), DIMENSION(jpi,jpj) :: utau_oce !: atm.-ocean surface i-stress (ocean referential) [N/m2] 45 REAL(wp), DIMENSION(jpi,jpj) :: vtau_oce !: atm.-ocean surface j-stress (ocean referential) [N/m2] 43 46 !! * Substitutions 44 47 # include "vectopt_loop_substitute.h90" 45 48 !!---------------------------------------------------------------------- 46 !! LIM 3.0, UCL-LOCEAN-IPSL (2006)47 !! $ 49 !! NEMO/LIM 3.0 , UCL-LOCEAN-IPSL (2008) 50 !! $Id: $ 48 51 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 49 52 !!---------------------------------------------------------------------- … … 51 54 CONTAINS 52 55 53 SUBROUTINE lim_sbc ( kt)56 SUBROUTINE lim_sbc_tau( kt, kcpl ) 54 57 !!------------------------------------------------------------------- 55 !! *** ROUTINE lim_sbc ***58 !! *** ROUTINE lim_sbc_tau *** 56 59 !! 57 !! ** Purpose : Update surface ocean boundary condition over areas 58 !! that are at least partially covered by sea-ice 60 !! ** Purpose : Update the ocean surface stresses due to the ice 59 61 !! 60 !! ** Action : - comput. of the momentum, heat and freshwater/salt 61 !! fluxes at the ice-ocean interface. 62 !! - Update 62 !! ** Action : - compute the ice-ocean stress depending on kcpl: 63 !! fluxes at the ice-ocean interface. 64 !! Case 0 : old LIM-3 way, computed at ice time-step only 65 !! Case 1 : at each ocean time step the stresses are computed 66 !! using the current ocean velocity (now) 67 !! Case 2 : combination of half case 0 + half case 1 68 !! 69 !! ** Outputs : - utau : sea surface i-stress (ocean referential) 70 !! - vtau : sea surface j-stress (ocean referential) 71 !! 72 !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 73 !! Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 74 !!--------------------------------------------------------------------- 75 INTEGER :: kt ! number of ocean iteration 76 INTEGER :: kcpl ! ice-ocean coupling indicator: =0 LIM-3 old case 77 ! ! =1 stresses computed using now ocean velocity 78 ! ! =2 combination of 0 and 1 cases 79 !! 80 INTEGER :: ji, jj ! dummy loop indices 81 REAL(wp) :: zfrldu, zfrldv ! lead fraction at U- & V-points 82 REAL(wp) :: ztglx , ztgly 83 REAL(wp) :: zat_u, zu_ico, zutaui, zu_u, zv_u, zmodu, zmod 84 REAL(wp) :: zat_v, zv_ico, zvtaui, zu_v, zv_v, zmodv, zsang 85 86 #if defined key_coupled 87 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb ! albedo of ice under overcast sky 88 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalbp ! albedo of ice under clear sky 89 #endif 90 REAL(wp), DIMENSION(jpi,jpj) :: ztio_u, ztio_v ! ocean stress below sea-ice 91 !!--------------------------------------------------------------------- 92 93 IF( kt == nit000 ) THEN 94 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM 3.0 sea-ice - surface ocean momentum fluxes' 96 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 97 ENDIF 98 99 SELECT CASE( kcpl ) 100 ! !--------------------------------! 101 CASE( 0 ) ! LIM 3 old stress computation ! (at ice timestep only) 102 ! !--------------------------------! 103 ! ... ice stress over ocean with a ice-ocean rotation angle 104 DO jj = 1, jpj 105 ! ... change the sinus angle sign in the south hemisphere 106 zsang = SIGN(1.e0, gphif(1,jj) ) * sangvg 107 DO ji = 1, jpi 108 ! ... ice velocity relative to the ocean 109 zu_ico = u_ice(ji,jj) - u_oce(ji,jj) 110 zv_ico = v_ice(ji,jj) - v_oce(ji,jj) 111 zmod = SQRT( zu_ico * zu_ico + zv_ico * zv_ico ) 112 ! quadratic drag formulation 113 ztglx = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico ) 114 ztgly = rhoco * zmod * ( cangvg * zv_ico + zsang * zu_ico ) 115 ! IMPORTANT 116 ! these lines are bound to prevent numerical oscillations 117 ! in the ice-ocean stress 118 ! They are physically ill-based. There is a cleaner solution 119 ! to try (remember discussion in Paris Gurvan) 120 ztio_u(ji,jj) = ztglx * exp( - zmod / 0.5 ) 121 ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 ) 122 ! 123 END DO 124 END DO 125 126 DO jj = 2, jpjm1 127 DO ji = fs_2, fs_jpim1 ! vertor opt. 128 ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 129 zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj ) ) 130 zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji ,jj+1) ) 131 ! update surface ocean stress 132 utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj) 133 vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj) 134 ! 135 END DO 136 END DO 137 138 ! boundary condition on the stress (utau,vtau) 139 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) 140 ! 141 ! !--------------------------------! 142 CASE( 1 ) ! Use the now velocity ! (at each ocean timestep) 143 ! !--------------------------------! 144 ! ... save the air-ocean stresses at ice time-step 145 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 146 utau_oce(:,:) = utau(:,:) 147 vtau_oce(:,:) = vtau(:,:) 148 ENDIF 149 ! ... ice stress over ocean with a ice-ocean rotation angle 150 DO jj = 2, jpjm1 151 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 152 DO ji = fs_2, fs_jpim1 153 ! computation of wind stress over ocean in X and Y direction 154 zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5 ! ice area at u and V-points 155 zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 156 157 zu_u = u_ice(ji,jj) - un(ji,jj,1) ! u ice-ocean velocity at U-point 158 zv_v = v_ice(ji,jj) - vn(ji,jj,1) ! v ice-ocean velocity at V-point 159 ! ! u ice-ocean velocity at V-point 160 zu_v = 0.25 * ( u_ice(ji,jj ) - un(ji,jj ,1) + u_ice(ji+1,jj ) - un(ji+1,jj ,1) & 161 & + u_ice(ji,jj-1) - un(ji,jj-1,1) + u_ice(ji+1,jj-1) - un(ji+1,jj-1,1) ) 162 ! ! v ice-ocean velocity at U-point 163 zv_u = 0.25 * ( v_ice(ji-1,jj+1) - vn(ji-1,jj+1,1) + v_ice(ji,jj+1) - vn(ji,jj+1,1) & 164 & + v_ice(ji-1,jj ) - vn(ji-1,jj ,1) + v_ice(ji,jj ) - vn(ji,jj ,1) ) 165 zmodu = SQRT( zu_u * zu_u + zv_u * zv_u ) ! module of the ice-ocean velocity at U-point 166 zmodv = SQRT( zu_v * zu_v + zv_v * zv_v ) ! module of the ice-ocean velocity at V-point 167 zutaui = rhoco * zmodu * ( cangvg * zu_u - zsang * zv_u ) 168 zvtaui = rhoco * zmodv * ( cangvg * zv_v + zsang * zu_v ) 169 170 utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui ! stress at the ocean surface 171 vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 172 END DO 173 END DO 174 ! boundary condition on the stress (utau,vtau) 175 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) 176 ! 177 ! !--------------------------------! 178 CASE( 2 ) ! mixed 0 and 2 cases ! (at each ocean timestep) 179 ! !--------------------------------! 180 ! ... save the air-ocean stresses at ice time-step 181 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 182 utau_oce(:,:) = utau(:,:) 183 vtau_oce(:,:) = vtau(:,:) 184 ENDIF 185 ! ... ice stress over ocean with a ice-ocean rotation angle 186 DO jj = 2, jpjm1 187 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 188 DO ji = fs_2, fs_jpim1 189 ! computation of wind stress over ocean in X and Y direction 190 zat_u = at_i(ji,jj) + at_i(ji+1,jj) * 0.5 ! ice area at u and V-points 191 zat_v = at_i(ji,jj) + at_i(ji,jj+1) * 0.5 192 193 !!gm bug mixing U and V points value below ====>>> to be corrected 194 zu_ico = u_ice(ji,jj) - 0.5 * ( un(ji,jj,1) - ssu_m(ji,jj) ) ! ice-oce velocity using un and ssu_m 195 zv_ico = v_ice(ji,jj) - 0.5 * ( vn(ji,jj,1) - ssu_m(ji,jj) ) 196 ! ! module of the ice-ocean velocity 197 zmod = SQRT( zu_ico * zu_ico + zv_ico * zv_ico ) 198 ! ! quadratic drag formulation with rotation 199 zutaui = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico ) 200 zvtaui = rhoco * zmod * ( cangvg * zv_ico + zsang * zu_ico ) 201 ! 202 utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui ! stress at the ocean surface 203 vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 204 END DO 205 END DO 206 ! boundary condition on the stress (utau,vtau) 207 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) 208 ! 209 END SELECT 210 IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, & 211 & tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask ) 212 ! 213 END SUBROUTINE lim_sbc_tau 214 215 216 SUBROUTINE lim_sbc_flx( kt ) 217 !!------------------------------------------------------------------- 218 !! *** ROUTINE lim_sbc_flx *** 219 !! 220 !! ** Purpose : Update the surface ocean boundary condition for heat 221 !! salt and mass over areas where sea-ice is non-zero 222 !! 223 !! ** Action : - computes the heat and freshwater/salt fluxes 224 !! at the ice-ocean interface. 225 !! - Update the ocean sbc 63 226 !! 64 227 !! ** Outputs : - qsr : sea heat flux: solar … … 66 229 !! - emp : freshwater budget: volume flux 67 230 !! - emps : freshwater budget: concentration/dillution 68 !! - utau : sea surface i-stress (ocean referential)69 !! - vtau : sea surface j-stress (ocean referential)70 231 !! 71 232 !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. … … 80 241 REAL(wp) :: zfons ! salt exchanges at the ice/ocean interface 81 242 REAL(wp) :: zpme ! freshwater exchanges at the ice/ocean interface 82 REAL(wp) :: zfrldu, zfrldv ! lead fraction at U- & V-points 83 REAL(wp) :: zu_io , zv_io ! 2 components of the ice-ocean velocity 84 REAL(wp) :: ztglx , ztgly ! temporary scalar 85 243 REAL(wp), DIMENSION(jpi,jpj) :: zfcm1 , zfcm2 ! solar/non solar heat fluxes 86 244 #if defined key_coupled 87 245 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb ! albedo of ice under overcast sky 88 246 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalbp ! albedo of ice under clear sky 89 247 #endif 90 REAL(wp) :: zsang, zmod91 REAL(wp), DIMENSION(jpi,jpj) :: ztio_u, ztio_v ! ocean stress below sea-ice92 REAL(wp), DIMENSION(jpi,jpj) :: zfcm1 , zfcm2 ! solar/non solar heat fluxes93 94 248 !!--------------------------------------------------------------------- 95 249 96 250 IF( kt == nit000 ) THEN 97 251 IF(lwp) WRITE(numout,*) 98 IF(lwp) WRITE(numout,*) 'lim_sbc : LIM 3.0 sea-ice - surface boundary condition'99 IF(lwp) WRITE(numout,*) '~~~~~~~ 252 IF(lwp) WRITE(numout,*) 'lim_sbc_flx : LIM 3.0 sea-ice - heat salt and mass ocean surface fluxes' 253 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 100 254 ENDIF 101 255 … … 169 323 ! ! fdtcn : turbulent oceanic heat flux 170 324 171 325 !!gm this IF prevents the vertorisation of the whole loop 172 326 IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN 173 327 WRITE(numout,*) ' lim_sbc : heat fluxes ' … … 198 352 WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 199 353 ENDIF 354 !!gm end 200 355 END DO 201 356 END DO … … 205 360 !------------------------------------------! 206 361 362 !!gm optimisation: this loop have to be merged with the previous one 207 363 DO jj = 1, jpj 208 364 DO ji = 1, jpi … … 254 410 END DO 255 411 256 emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 257 IF (num_sal.eq.2) THEN 258 !In case of variable salinity the salt flux has to be accounted for differently 259 ! Brine drainage has to be added 412 IF( num_sal == 2 ) THEN ! variable ice salinity: brine drainage included in the salt flux 260 413 emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 414 ELSE ! constant ice salinity: 415 emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 261 416 ENDIF 262 417 263 418 IF( lk_dynspg_rl ) emp (:,:) = emps(:,:) ! rigid-lid formulation : emp = emps 264 265 !------------------------------------------!266 ! momentum flux at the ocean surface !267 !------------------------------------------!268 269 IF ( ln_limdyn ) THEN ! Update the stress over ice-over area (only in ice-dynamic case)270 ! ! otherwise the atmosphere-ocean stress is used everywhere271 272 ! ... ice stress over ocean with a ice-ocean rotation angle273 DO jj = 1, jpj274 ! ... change the sinus angle sign in the south hemisphere275 zsang = SIGN(1.e0, gphif(1,jj) ) * sangvg276 DO ji = 1, jpi277 ! ... ice velocity relative to the ocean278 zu_io = u_ice(ji,jj) - u_oce(ji,jj)279 zv_io = v_ice(ji,jj) - v_oce(ji,jj)280 zmod = SQRT( zu_io * zu_io + zv_io * zv_io )281 ! quadratic drag formulation282 ztglx = rhoco * zmod * ( cangvg * zu_io - zsang * zv_io )283 ztgly = rhoco * zmod * ( cangvg * zv_io + zsang * zu_io )284 ! IMPORTANT285 ! these lines are bound to prevent numerical oscillations286 ! in the ice-ocean stress287 ! They are physically ill-based. There is a cleaner solution288 ! to try (remember discussion in Paris Gurvan)289 ztio_u(ji,jj) = ztglx * exp( - zmod / 0.5 )290 ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 )291 !292 END DO293 END DO294 295 DO jj = 2, jpjm1296 DO ji = fs_2, fs_jpim1 ! vertor opt.297 ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)298 zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj ) )299 zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji ,jj+1) )300 ! update surface ocean stress301 utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj)302 vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj)303 !304 END DO305 END DO306 307 ! boundary condition on the stress (utau,vtau)308 CALL lbc_lnk( utau, 'U', -1. )309 CALL lbc_lnk( vtau, 'V', -1. )310 311 ENDIF312 419 313 420 !-----------------------------------------------! … … 331 438 332 439 IF(ln_ctl) THEN 333 CALL prt_ctl(tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ') 334 CALL prt_ctl(tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=emps , clinfo2=' emps : ') 335 CALL prt_ctl(tab2d_1=utau , clinfo1=' lim_sbc: utau : ', mask1=umask, & 336 & tab2d_2=vtau , clinfo2=' vtau : ' , mask2=vmask ) 337 CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ') 338 CALL prt_ctl(tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl) 440 CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ' ) 441 CALL prt_ctl( tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=emps, clinfo2=' emps : ' ) 442 CALL prt_ctl( tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ' ) 443 CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl ) 339 444 ENDIF 340 341 END SUBROUTINE lim_sbc 445 ! 446 END SUBROUTINE lim_sbc_flx 447 342 448 343 449 #else
Note: See TracChangeset
for help on using the changeset viewer.