Changeset 716 for trunk/NEMO
- Timestamp:
- 2007-10-11T17:13:00+02:00 (17 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC/limdyn.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC/limdyn.F90
r699 r716 34 34 REAL(wp) :: rone = 1.e0 ! constant value 35 35 36 37 # include "vectopt_loop_substitute.h90" 38 36 39 !!---------------------------------------------------------------------- 37 40 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) … … 64 67 REAL(wp) :: & 65 68 ztairx, ztairy, & ! tempory scalars 66 zsang , z mod, &69 zsang , zrhomod, & 67 70 ztglx , ztgly , & 68 71 zt11, zt12, zt21, zt22 , & 69 72 zustm, zsfrld, zsfrldm4, & 70 73 zu_ice, zv_ice, ztair2 74 REAL(wp),DIMENSION(jpi,jpj) :: zmod 71 75 REAL(wp),DIMENSION(jpj) :: & 72 76 zind, & ! i-averaged indicator of sea-ice … … 82 86 hicm(:,:) = ( 1.0 - frld(:,:) ) * hicif(:,:) 83 87 84 u_ oce(:,:) = u_io(:,:) * tmu(:,:)85 v_ oce(:,:) = v_io(:,:) * tmu(:,:)88 u_io(:,:) = u_io(:,:) * tmu(:,:) 89 v_io(:,:) = v_io(:,:) * tmu(:,:) 86 90 87 91 ! ! Rheology (ice dynamics) … … 156 160 157 161 IF(ln_ctl) THEN 158 CALL prt_ctl(tab2d_1=u_ oce , clinfo1=' lim_dyn : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce:')159 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice:')162 CALL prt_ctl(tab2d_1=u_io , clinfo1=' lim_dyn : u_io :', tab2d_2=v_io , clinfo2=' v_io :') 163 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_dyn : u_ice:', tab2d_2=v_ice, clinfo2=' v_ice:') 160 164 ENDIF 161 165 162 166 ! ! Ice-Ocean stress 163 ! ! ================ 164 DO jj = 2, jpjm1 165 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 166 DO ji = 2, jpim1 167 ! computation of wind stress over ocean in X and Y direction 168 #if defined key_coupled && defined key_lim_cp1 169 ztairx = frld(ji-1,jj ) * gtaux(ji-1,jj ) + frld(ji,jj ) * gtaux(ji,jj ) & 170 & + frld(ji-1,jj-1) * gtaux(ji-1,jj-1) + frld(ji,jj-1) * gtaux(ji,jj-1) 171 172 ztairy = frld(ji-1,jj ) * gtauy(ji-1,jj ) + frld(ji,jj ) * gtauy(ji,jj ) & 173 & + frld(ji-1,jj-1) * gtauy(ji-1,jj-1) + frld(ji,jj-1) * gtauy(ji,jj-1) 174 #else 175 zsfrld = frld(ji,jj) + frld(ji-1,jj) + frld(ji-1,jj-1) + frld(ji,jj-1) 176 ztairx = zsfrld * gtaux(ji,jj) 177 ztairy = zsfrld * gtauy(ji,jj) 178 #endif 179 zsfrldm4 = 4 - frld(ji,jj) - frld(ji-1,jj) - frld(ji-1,jj-1) - frld(ji,jj-1) 180 zu_ice = u_ice(ji,jj) - u_oce(ji,jj) 181 zv_ice = v_ice(ji,jj) - v_oce(ji,jj) 182 zmod = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 183 ztglx = zsfrldm4 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 184 ztgly = zsfrldm4 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 185 186 tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 4 * rau0 ) 187 tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 4 * rau0 ) 167 ! ! ================ 168 DO jj = 1, jpj 169 DO ji = 1, jpi 170 !! zsang = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg ! do the full loop and avoid lbc_lnk 171 zsang = SIGN(1.e0, gphif(ji,jj) ) * sangvg 172 zu_ice = u_ice(ji,jj) - u_io(ji,jj) 173 zv_ice = v_ice(ji,jj) - v_io(ji,jj) 174 zrhomod = zu_ice * zu_ice + zv_ice * zv_ice 175 zmod (ji,jj) = zrhomod 176 zrhomod = rhoco * SQRT( zrhomod ) 177 ftaux(ji,jj) = zrhomod * ( cangvg * zu_ice - zsang * zv_ice ) 178 ftauy(ji,jj) = zrhomod * ( cangvg * zv_ice + zsang * zu_ice ) 188 179 END DO 189 180 END DO 190 181 191 182 ! computation of friction velocity 192 183 DO jj = 2, jpjm1 193 DO ji = 2, jpim1 194 195 zu_ice = u_ice(ji-1,jj-1) - u_oce(ji-1,jj-1) 196 zv_ice = v_ice(ji-1,jj-1) - v_oce(ji-1,jj-1) 197 zt11 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 198 199 zu_ice = u_ice(ji-1,jj) - u_oce(ji-1,jj) 200 zv_ice = v_ice(ji-1,jj) - v_oce(ji-1,jj) 201 zt12 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 202 203 zu_ice = u_ice(ji,jj-1) - u_oce(ji,jj-1) 204 zv_ice = v_ice(ji,jj-1) - v_oce(ji,jj-1) 205 zt21 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 206 207 zu_ice = u_ice(ji,jj) - u_oce(ji,jj) 208 zv_ice = v_ice(ji,jj) - v_oce(ji,jj) 209 zt22 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 210 211 ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj) 212 213 zustm = ( 1 - frld(ji,jj) ) * 0.25 * ( zt11 + zt12 + zt21 + zt22 ) & 214 & + frld(ji,jj) * SQRT( ztair2 ) 215 216 ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj) 184 DO ji = fs_2, fs_jpim1 ! vector opt. 185 ust2s(ji,jj) = 0.25 * cw * ( zmod(ji,jj+1) + zmod(ji+1,jj+1) + & 186 & zmod(ji,jj ) + zmod(ji+1,jj ) ) * tms(ji,jj) 217 187 END DO 218 188 END DO … … 220 190 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 221 191 192 ftaux(:,:) = gtaux(:,:) 193 ftauy(:,:) = gtauy(:,:) 194 222 195 DO jj = 2, jpjm1 223 196 DO ji = 2, jpim1 224 #if defined key_coupled && defined key_lim_cp1 225 tio_u(ji,jj) = - ( gtaux(ji ,jj ) + gtaux(ji-1,jj ) & 226 & + gtaux(ji-1,jj-1) + gtaux(ji ,jj-1) ) / ( 4 * rau0 ) 227 228 tio_v(ji,jj) = - ( gtauy(ji ,jj ) + gtauy(ji-1,jj ) & 229 & + gtauy(ji-1,jj-1) + gtauy(ji ,jj-1) ) / ( 4 * rau0 ) 230 #else 231 tio_u(ji,jj) = - gtaux(ji,jj) / rau0 232 tio_v(ji,jj) = - gtauy(ji,jj) / rau0 233 #endif 234 ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj) 235 zustm = SQRT( ztair2 ) 236 237 ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj) 197 ztair2 = gtaux(ji ,jj+1) * gtaux(ji ,jj+1) + gtauy(ji ,jj+1) * gtauy(ji ,jj+1) & 198 & + gtaux(ji+1,jj+1) * gtaux(ji+1,jj+1) + gtauy(ji+1,jj+1) * gtauy(ji+1,jj+1) & 199 & + gtaux(ji ,jj ) * gtaux(ji ,jj ) + gtauy(ji ,jj ) * gtauy(ji ,jj ) & 200 & + gtaux(ji+1,jj ) * gtaux(ji+1,jj ) + gtauy(ji+1,jj ) * gtauy(ji+1,jj ) 201 202 ust2s(ji,jj) = 0.25 / rau0 * SQRT( ztair2 ) * tms(ji,jj) 238 203 END DO 239 204 END DO … … 242 207 243 208 CALL lbc_lnk( ust2s, 'T', 1. ) ! T-point 244 CALL lbc_lnk( tio_u, 'I', -1. ) ! I-point (i.e. ice U-V point)245 CALL lbc_lnk( tio_v, 'I', -1. ) ! I-point (i.e. ice U-V point)246 209 247 210 IF(ln_ctl) THEN 248 CALL prt_ctl(tab2d_1= tio_u , clinfo1=' lim_dyn : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v:')211 CALL prt_ctl(tab2d_1=ftaux , clinfo1=' lim_dyn : ftaux :', tab2d_2=ftauy , clinfo2=' ftauy :') 249 212 CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn : ust2s :') 250 213 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
