Changeset 1526
- Timestamp:
- 2009-07-22T15:51:02+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limsbc.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limsbc.F90
r1525 r1526 40 40 REAL(wp) :: rone = 1.e0 41 41 42 REAL(wp), DIMENSION(jpi,jpj) :: utau_oce, vtau_oce !: air-ocean surface i- & j-stress [N/m2]42 REAL(wp), DIMENSION(jpi,jpj) :: utau_oce, vtau_oce !: air-ocean surface i- & j-stress [N/m2] 43 43 REAL(wp), DIMENSION(jpi,jpj) :: tmod_io !: modulus of the ice-ocean relative velocity [m/s] 44 REAL(wp), DIMENSION(jpi,jpj) :: ssu_mb , ssv_mb !: before mean ocean surface currents [m/s] 44 REAL(wp), DIMENSION(jpi,jpj) :: ssu_mb , ssv_mb !: before mean ocean surface currents [m/s] 45 45 46 !! * Substitutions 46 47 # include "vectopt_loop_substitute.h90" 47 48 !!---------------------------------------------------------------------- 48 !! NEMO/LIM 3. 0 , UCL-LOCEAN-IPSL (2008)49 !! NEMO/LIM 3.2 , UCL-LOCEAN-IPSL (2009) 49 50 !! $Id$ 50 51 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 77 78 ! ! =2 combination of 0 and 1 cases 78 79 !! 79 INTEGER :: ji, jj ! dummy loop indices 80 REAL(wp) :: zfrldu, zfrldv ! lead fraction at U- & V-points 81 REAL(wp) :: zat_u, zu_ico, zutaui, zu_u 82 REAL(wp) :: zat_v, zv_ico, zvtaui, zv_v, zsang 83 REAL(wp) :: zu_ij, zu_im1j, zv_ij, zv_ijm1 84 80 INTEGER :: ji, jj ! dummy loop indices 81 REAL(wp) :: zfrldu, zat_u, zu_ico, zutaui, zu_u, zu_ij, zu_im1j ! temporary scalar 82 REAL(wp) :: zfrldv, zat_v, zv_ico, zvtaui, zv_v, zv_ij, zv_ijm1 ! - - 83 REAL(wp) :: zsang ! - - 84 REAL(wp), DIMENSION(jpi,jpj) :: ztio_u, ztio_v ! ocean stress below sea-ice 85 85 #if defined key_coupled 86 86 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb ! albedo of ice under overcast sky 87 87 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalbp ! albedo of ice under clear sky 88 88 #endif 89 REAL(wp), DIMENSION(jpi,jpj) :: ztio_u, ztio_v ! ocean stress below sea-ice 90 !!--------------------------------------------------------------------- 89 !!--------------------------------------------------------------------- 91 90 92 91 IF( kt == nit000 ) THEN … … 100 99 CASE( 0 ) ! LIM 3 old stress computation ! (at ice timestep only) 101 100 ! !--------------------------------! 102 DO jj = 2, jpjm1 ! ...modulus of the ice-ocean velocity101 DO jj = 2, jpjm1 !* modulus of the ice-ocean velocity 103 102 DO ji = fs_2, fs_jpim1 104 103 zu_ij = u_ice(ji ,jj) - ssu_m(ji ,jj) ! (i ,j) … … 110 109 END DO 111 110 END DO 112 CALL lbc_lnk( tmod_io, 'T', 1. ) 113 ! ... ice stress over ocean with a ice-ocean rotation angle 111 CALL lbc_lnk( tmod_io, 'T', 1. ) ! lateral boundary condition 112 ! 113 ! !* ice stress over ocean with a ice-ocean rotation angle 114 114 DO jj = 1, jpjm1 115 ! ... change the sinus angle sign in the south hemisphere 116 zsang = SIGN(1.e0, gphif(1,jj) ) * sangvg 115 zsang = SIGN( 1.e0, gphif(1,jj) ) * sangvg ! change the sinus angle sign in the south hemisphere 117 116 DO ji = 1, fs_jpim1 118 ! ... ice velocity relative to the ocean 119 zu_u = u_ice(ji,jj) - u_oce(ji,jj) 117 zu_u = u_ice(ji,jj) - u_oce(ji,jj) ! ice velocity relative to the ocean 120 118 zv_v = v_ice(ji,jj) - v_oce(ji,jj) 121 ! quadratic drag formulation119 ! ! quadratic drag formulation with rotation 122 120 !!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 123 121 zutaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v ) 124 122 zvtaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u ) 125 ! IMPORTANT 126 ! these lines are bound to prevent numerical oscillations 127 ! in the ice-ocean stress 128 ! They are physically ill-based. There is a cleaner solution 129 ! to try (remember discussion in Paris Gurvan) 130 ztio_u(ji,jj) = zutaui * exp( - ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) ) 131 ztio_v(ji,jj) = zvtaui * exp( - ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) ) 132 ! 133 END DO 134 END DO 135 136 DO jj = 2, jpjm1 123 ! ! bound for too large stress values 124 ! IMPORTANT: the exponential below prevents numerical oscillations in the ice-ocean stress 125 ! This is not physically based. A cleaner solution is offer in CASE kcpl=2 126 ztio_u(ji,jj) = zutaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) ) 127 ztio_v(ji,jj) = zvtaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) ) 128 END DO 129 END DO 130 DO jj = 2, jpjm1 ! stress at the surface of the ocean 137 131 DO ji = fs_2, fs_jpim1 ! vertor opt. 138 ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 139 zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj ) ) 132 zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj ) ) ! open-ocean fraction at U- & V-points (from T-point values) 140 133 zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji ,jj+1) ) 141 ! update surface ocean stress134 ! ! update surface ocean stress 142 135 utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj) 143 136 vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj) 144 ! 145 END DO 146 END DO 147 148 ! boundary condition on the stress (utau,vtau) 149 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) 137 END DO 138 END DO 139 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition 140 150 141 ! 151 142 ! !--------------------------------! … … 153 144 ! !--------------------------------! 154 145 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 155 utau_oce(:,:) = utau(:,:) ! ...save the air-ocean stresses at ice time-step146 utau_oce(:,:) = utau(:,:) !* save the air-ocean stresses at ice time-step 156 147 vtau_oce(:,:) = vtau(:,:) 157 DO jj = 2, jpjm1 ! ...modulus of the ice-ocean velocity148 DO jj = 2, jpjm1 !* modulus of the ice-ocean velocity 158 149 DO ji = fs_2, fs_jpim1 159 150 zu_ij = u_ice(ji ,jj) - ssu_m(ji ,jj) ! (i ,j) … … 165 156 END DO 166 157 END DO 167 CALL lbc_lnk( tmod_io, 'T', 1. )158 CALL lbc_lnk( tmod_io, 'T', 1. ) ! lateral boundary condition 168 159 ENDIF 169 ! ... ice stress over ocean with a ice-ocean rotation angle170 DO jj = 2, jpjm1171 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 160 ! 161 DO jj = 2, jpjm1 !* ice stress over ocean with a ice-ocean rotation angle 162 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg ! change the sinus angle sign in the south hemisphere 172 163 DO ji = fs_2, fs_jpim1 173 ! computation of wind stress over ocean in X and Y direction 174 zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5 ! ice area at u and V-points 164 zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5 ! ice area at u and V-points 175 165 zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 176 177 zu_u = u_ice(ji,jj) - ub(ji,jj,1) ! u ice-ocean velocity at U-point178 zv_v = v_ice(ji,jj) - vb(ji,jj,1) ! v ice-ocean velocity at V-point179 ! 166 ! ! (u,v) ice-ocean velocity at (U,V)-point, resp. 167 zu_u = u_ice(ji,jj) - ub(ji,jj,1) 168 zv_v = v_ice(ji,jj) - vb(ji,jj,1) 169 ! ! quadratic drag formulation with rotation 180 170 !!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 181 171 zutaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v ) 182 172 zvtaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u ) 183 184 utau(ji,jj) = ( 1.- zat_u ) * utau_oce(ji,jj) + zat_u * zutaui ! stress at the ocean surface173 ! ! stress at the ocean surface 174 utau(ji,jj) = ( 1.- zat_u ) * utau_oce(ji,jj) + zat_u * zutaui 185 175 vtau(ji,jj) = ( 1.- zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 186 ! ! u ice-ocean velocity at V-point 187 END DO 188 END DO 189 ! boundary condition on the stress (utau,vtau) 190 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) 176 END DO 177 END DO 178 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition 179 191 180 ! 192 181 ! !--------------------------------! … … 194 183 ! !--------------------------------! 195 184 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 196 utau_oce(:,:) = utau (:,:) ! ...save the air-ocean stresses at ice time-step185 utau_oce(:,:) = utau (:,:) !* save the air-ocean stresses at ice time-step 197 186 vtau_oce(:,:) = vtau (:,:) 198 ssu_mb (:,:) = ssu_m(:,:) ! ...save the ice-ocean velocity at ice time-step187 ssu_mb (:,:) = ssu_m(:,:) !* save the ice-ocean velocity at ice time-step 199 188 ssv_mb (:,:) = ssv_m(:,:) 200 DO jj = 2, jpjm1 ! ...modulus of the ice-ocean velocity189 DO jj = 2, jpjm1 !* modulus of the ice-ocean velocity 201 190 DO ji = fs_2, fs_jpim1 202 191 zu_ij = u_ice(ji ,jj) - ssu_m(ji ,jj) ! (i ,j) … … 208 197 END DO 209 198 END DO 210 CALL lbc_lnk( tmod_io, 'T', 1. )199 CALL lbc_lnk( tmod_io, 'T', 1. ) 211 200 ENDIF 212 ! ... ice stress over ocean with a ice-ocean rotation angle 213 DO jj = 2, jpjm1 201 DO jj = 2, jpjm1 !* ice stress over ocean with a ice-ocean rotation angle 214 202 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 215 203 DO ji = fs_2, fs_jpim1 216 ! computation of wind stress over ocean in X and Y direction217 204 zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5 ! ice area at u and V-points 218 205 zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 219 220 !!gm bug mixing U and V points value below ====>>> to be corrected 221 zu_ico = u_ice(ji,jj) - 0.5 * ( ub(ji,jj,1) - ssu_mb(ji,jj) ) ! ice-oce velocity using un and ssu_mb 206 ! 207 zu_ico = u_ice(ji,jj) - 0.5 * ( ub(ji,jj,1) - ssu_mb(ji,jj) ) ! ice-oce velocity using ub and ssu_mb 222 208 zv_ico = v_ice(ji,jj) - 0.5 * ( vb(ji,jj,1) - ssv_mb(ji,jj) ) 223 209 ! ! quadratic drag formulation with rotation … … 230 216 END DO 231 217 END DO 232 ! boundary condition on the stress (utau,vtau) 233 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) 218 CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition 234 219 ! 235 220 END SELECT 221 236 222 IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, & 237 223 & tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask ) … … 389 375 !------------------------------------------! 390 376 391 377 !!gm optimisation: this loop have to be merged with the previous one 392 378 DO jj = 1, jpj 393 379 DO ji = 1, jpi … … 475 461 END SUBROUTINE lim_sbc_flx 476 462 477 478 463 #else 479 464 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
