Changeset 1530 for trunk/NEMO/LIM_SRC_2
- Timestamp:
- 2009-07-23T18:22:51+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_2/limadv_2.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_2/limadv_2.F90
r1529 r1530 4 4 !! LIM 2.0 sea-ice model : sea-ice advection 5 5 !!====================================================================== 6 !! History : OPA ! 2000-01 (LIM) Original code 7 !! ! 2001-05 (G. Madec, R. Hordoir) Doctor norm 8 !! NEMO 1.0 ! 2003-10 (C. Ethe) F90, module 9 !! - ! 2003-12 (R. Hordoir, G. Madec) mpp 10 !! 3.2 ! 2009-06 (F. Dupont) correct a error in the North fold b. c. 11 !!-------------------------------------------------------------------- 6 12 #if defined key_lim2 7 13 !!---------------------------------------------------------------------- … … 11 17 !! lim_adv_y_2 : advection of sea ice on y axis 12 18 !!---------------------------------------------------------------------- 13 !! * Modules used14 19 USE dom_oce 15 20 USE dom_ice_2 16 USE in_out_manager ! I/O manager17 21 USE ice_2 18 22 USE lbclnk 19 USE prtctl ! Print control 23 USE in_out_manager ! I/O manager 24 USE prtctl ! Print control 20 25 21 26 IMPLICIT NONE 22 27 PRIVATE 23 28 24 !! * Routine accessibility 25 PUBLIC lim_adv_x_2 ! called by lim_trp 26 PUBLIC lim_adv_y_2 ! called by lim_trp 27 29 PUBLIC lim_adv_x_2 ! called by lim_trp 30 PUBLIC lim_adv_y_2 ! called by lim_trp 31 32 REAL(wp) :: epsi20 = 1.e-20 ! constant values 33 REAL(wp) :: rzero = 0.e0 ! - - 34 REAL(wp) :: rone = 1.e0 ! - - 35 28 36 !! * Substitutions 29 37 # include "vectopt_loop_substitute.h90" 30 31 !! * Module variables 32 REAL(wp) :: & ! constant values 33 epsi20 = 1e-20 , & 34 rzero = 0.e0 , & 35 rone = 1.e0 36 !!---------------------------------------------------------------------- 37 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 38 !!---------------------------------------------------------------------- 39 !! NEMO/LIM 3.2, UCL-LOCEAN-IPSL (2009) 38 40 !! $Id$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 42 !!---------------------------------------------------------------------- 41 43 … … 48 50 !! 49 51 !! ** purpose : Computes and adds the advection trend to sea-ice 50 !! variable on xaxis52 !! variable on i-axis 51 53 !! 52 !! ** method : Uses Prather second order scheme that advects 53 !! tracers but also theirquadratic forms. The method preserves 54 !! tracer structures by conserving second order moments. 55 !! Ref.: "Numerical Advection by Conservation of Second Order 56 !! Moments", JGR, VOL. 91. NO. D6. PAGES 6671-6681. MAY 20, 1986 57 !! 58 !! History : 59 !! ! 00-01 (LIM) 60 !! ! 01-05 (G. Madec, R. Hordoir) opa norm 61 !! ! 03-10 (C. Ethe) F90, module 62 !! ! 03-12 (R. Hordoir, G. Madec) mpp 54 !! ** method : Uses Prather second order scheme that advects tracers 55 !! but also theirquadratic forms. The method preserves 56 !! tracer structures by conserving second order moments. 57 !! 58 !! Reference: Prather, 1986, JGR, 91, D6. 6671-6681. 63 59 !!-------------------------------------------------------------------- 64 !! * Arguments 65 REAL(wp) , INTENT(in) :: & 66 pdf , & ! ??? 67 pcrh ! = 1. : lim_adv_x is called before lim_adv_y 68 ! ! = 0. : lim_adv_x is called after lim_adv_y 69 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: & 70 put ! i-direction ice velocity at ocean U-point (m/s) 71 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & 72 ps0 , psm , & ! ??? 73 psx , psy , & ! ??? 74 psxx, psyy, psxy 75 76 !! * Local declarations 77 INTEGER :: ji, jj ! dummy loop indices 78 REAL(wp) :: & 79 zrdt, zslpmax, ztemp, zin0, & ! temporary scalars 80 zs1max, zs1new, zs2new, & ! " " 81 zalf, zalfq, zalf1, zalf1q, & ! " " 82 zbt , zbt1 ! " " 83 REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace 84 zf0 , zfx , zfy , zbet, & ! " " 85 zfxx, zfyy, zfxy, & ! " " 86 zfm, zalg, zalg1, zalg1q ! " " 60 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for the time step 61 REAL(wp) , INTENT(in ) :: pcrh ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 62 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 63 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psm ! area 64 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ps0 ! field to be advected 65 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments 66 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 67 !! 68 INTEGER :: ji, jj ! dummy loop indices 69 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp, zin0 ! temporary scalars 70 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 71 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 72 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 73 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - 74 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 87 75 !--------------------------------------------------------------------- 88 76 89 77 ! Limitation of moments. 90 78 91 zrdt = rdt_ice * pdf! If ice drift field is too fast, use an appropriate time step for advection.79 zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 92 80 93 81 DO jj = 1, jpj … … 99 87 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) ) ) 100 88 zin0 = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask 101 89 ! 102 90 ps0 (ji,jj) = zslpmax 103 91 psx (ji,jj) = zs1new * zin0 … … 120 108 zalf1 = 1.0 - zalf 121 109 zalf1q = zalf1 * zalf1 110 ! 122 111 zfm (ji,jj) = zalf * psm(ji,jj) 123 112 zf0 (ji,jj) = zalf * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) ) ) … … 127 116 zfxy(ji,jj) = zalfq * psxy(ji,jj) 128 117 zfyy(ji,jj) = zalf * psyy(ji,jj) 129 118 ! 130 119 ! Readjust moments remaining in the box. 131 120 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj) … … 181 170 zalf1 = 1.0 - zalf 182 171 ztemp = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj) 183 ps0(ji,jj) = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj) 184 psx(ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 185 psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj) & 172 ! 173 ps0 (ji,jj) = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj) 174 psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 175 psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj) & 186 176 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 187 177 & + zbt1 * psxx(ji,jj) … … 202 192 zalf1 = 1.0 - zalf 203 193 ztemp = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 204 ps0(ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 205 psx(ji,jj) = zbt * psx(ji,jj) & 206 & + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 207 psxx(ji,jj) = zbt * psxx(ji,jj) & 208 & + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj) & 209 & + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 210 psxy(ji,jj) = zbt * psxy(ji,jj) & 211 & + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) & 212 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) ) ) 213 psy(ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 214 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) 194 ! 195 ps0(ji,jj) = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 196 psx(ji,jj) = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 197 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj) & 198 & + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) & 199 & + ( zalf1 - zalf ) * ztemp ) ) 200 psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) & 201 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) ) ) 202 psy(ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 203 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) 215 204 END DO 216 205 END DO 217 206 218 207 !-- Lateral boundary conditions 219 CALL lbc_lnk( psm , 'T', 1. ) 220 CALL lbc_lnk( ps0 , 'T', 1. ) 221 CALL lbc_lnk( psx , 'T', -1. ) 222 CALL lbc_lnk( psxx, 'T', 1. ) 223 CALL lbc_lnk( psy , 'T', -1. ) 224 CALL lbc_lnk( psyy, 'T', 1. ) 208 CALL lbc_lnk( psm , 'T', 1. ) ; CALL lbc_lnk( ps0 , 'T', 1. ) 209 CALL lbc_lnk( psx , 'T', -1. ) ; CALL lbc_lnk( psy , 'T', -1. ) ! caution gradient ==> the sign changes 210 CALL lbc_lnk( psxx, 'T', 1. ) ; CALL lbc_lnk( psyy, 'T', 1. ) 225 211 CALL lbc_lnk( psxy, 'T', 1. ) 226 212 … … 231 217 CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :') 232 218 ENDIF 233 219 ! 234 220 END SUBROUTINE lim_adv_x_2 235 221 … … 241 227 !! 242 228 !! ** purpose : Computes and adds the advection trend to sea-ice 243 !! variable on yaxis229 !! variable on j-axis 244 230 !! 245 231 !! ** method : Uses Prather second order scheme that advects tracers 246 !! but also their quadratic forms. The method preserves tracer247 !! structures by conserving second order moments.232 !! but also their quadratic forms. The method preserves 233 !! tracer structures by conserving second order moments. 248 234 !! 249 !! History : 250 !! 1.0 ! 00-01 (LIM) 251 !! ! 01-05 (G. Madec, R. Hordoir) opa norm 252 !! 2.0 ! 03-10 (C. Ethe) F90, module 253 !! ! 03-12 (R. Hordoir, G. Madec) mpp 235 !! Reference: Prather, 1986, JGR, 91, D6. 6671-6681. 254 236 !!--------------------------------------------------------------------- 255 !! * Arguments 256 REAL(wp), INTENT(in) :: & 257 pdf, & ! ??? 258 pcrh ! = 1. : lim_adv_x is called before lim_adv_y 259 ! ! = 0. : lim_adv_x is called after lim_adv_y 260 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: & 261 pvt ! j-direction ice velocity at ocean V-point (m/s) 262 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & 263 psm , ps0 , psx , psy, & 264 psxx, psyy, psxy 265 266 !! * Local Variables 267 INTEGER :: ji, jj ! dummy loop indices 268 REAL(wp) :: & 269 zrdt, zslpmax, zin0, ztemp, & ! temporary scalars 270 zs1max, zs1new, zs2new, & ! " " 271 zalf, zalfq, zalf1, zalf1q, & ! " " 272 zbt , zbt1 ! 273 REAL(wp), DIMENSION(jpi,jpj) :: & 274 zf0 , zfx , zfy , & ! temporary workspace 275 zfxx, zfyy, zfxy, & ! " " 276 zfm , zbet, & ! " " 277 zalg, zalg1, zalg1q ! " " 237 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for the time step 238 REAL(wp) , INTENT(in ) :: pcrh ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 239 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 240 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psm ! area 241 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ps0 ! field to be advected 242 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments 243 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 244 !! 245 INTEGER :: ji, jj ! dummy loop indices 246 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp, zin0 ! temporary scalars 247 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 248 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 249 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 250 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - 251 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 278 252 !--------------------------------------------------------------------- 279 253 … … 282 256 zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 283 257 284 285 286 287 288 289 290 291 292 ps0 (ji,jj) = zslpmax293 psx (ji,jj) = psx (ji,jj) * zin0294 psxx(ji,jj) = psxx(ji,jj) * zin0295 psy (ji,jj) = zs1new* zin0296 psyy(ji,jj) = zs2new * zin0297 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) )* zin0298 END DO299 END DO300 301 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 302 psm (:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20)303 304 ! Calculate fluxes and moments between boxes j<-->j+1 305 DO jj = 1, jpj306 DO ji = 1, jpi307 ! Flux from j to j+1 WHEN v GT 0308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 258 DO jj = 1, jpj 259 DO ji = 1, jpi 260 zslpmax = MAX( rzero, ps0(ji,jj) ) 261 zs1max = 1.5 * zslpmax 262 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 263 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 264 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) ) ) 265 zin0 = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask 266 ! 267 ps0 (ji,jj) = zslpmax 268 psx (ji,jj) = psx (ji,jj) * zin0 269 psxx(ji,jj) = psxx(ji,jj) * zin0 270 psy (ji,jj) = zs1new * zin0 271 psyy(ji,jj) = zs2new * zin0 272 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0 273 END DO 274 END DO 275 276 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 277 psm(:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 278 279 ! Calculate fluxes and moments between boxes j<-->j+1 280 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 281 DO ji = 1, jpi 282 zbet(ji,jj) = MAX( rzero, SIGN( rone, pvt(ji,jj) ) ) 283 zalf = MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 284 zalfq = zalf * zalf 285 zalf1 = 1.0 - zalf 286 zalf1q = zalf1 * zalf1 287 zfm (ji,jj) = zalf * psm(ji,jj) 288 zf0 (ji,jj) = zalf * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj) + (zalf1-zalf) * psyy(ji,jj) ) ) 289 zfy (ji,jj) = zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 290 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj) 291 zfx (ji,jj) = zalf * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 292 zfxy(ji,jj) = zalfq * psxy(ji,jj) 293 zfxx(ji,jj) = zalf * psxx(ji,jj) 294 ! 295 ! Readjust moments remaining in the box. 296 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj) 297 ps0 (ji,jj) = ps0 (ji,jj) - zf0(ji,jj) 298 psy (ji,jj) = zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 299 psyy(ji,jj) = zalf1 * zalf1q * psyy(ji,jj) 300 psx (ji,jj) = psx (ji,jj) - zfx(ji,jj) 301 psxx(ji,jj) = psxx(ji,jj) - zfxx(ji,jj) 302 psxy(ji,jj) = zalf1q * psxy(ji,jj) 303 END DO 304 END DO 305 ! 306 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 307 DO ji = 1, jpi 308 zalf = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 309 zalg (ji,jj) = zalf 310 zalfq = zalf * zalf 311 zalf1 = 1.0 - zalf 312 zalg1 (ji,jj) = zalf1 313 zalf1q = zalf1 * zalf1 314 zalg1q(ji,jj) = zalf1q 315 zfm (ji,jj) = zfm (ji,jj) + zalf * psm(ji,jj+1) 316 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 317 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 318 zfyy (ji,jj) = zfyy(ji,jj) + zalf * zalfq * psyy(ji,jj+1) 319 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 320 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 321 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1) 322 END DO 323 END DO 350 324 351 ! Readjust moments remaining in the box. 352 DO jj = 2, jpj 353 DO ji = 1, jpi 354 zbt = zbet(ji,jj-1) 355 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 356 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 357 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 358 psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 359 psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj) 360 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 361 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 362 psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 363 END DO 364 END DO 365 366 ! Put the temporary moments into appropriate neighboring boxes. 367 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 368 DO ji = 1, jpi 369 zbt = zbet(ji,jj-1) 370 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 371 psm(ji,jj) = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 372 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj) 373 zalf1 = 1.0 - zalf 374 ztemp = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 375 ps0(ji,jj) = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj) 376 377 psy(ji,jj) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) & 378 & + zbt1 * psy(ji,jj) 379 380 psyy(ji,jj) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj) & 381 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 382 & + zbt1 * psyy(ji,jj) 383 384 psxy(ji,jj) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj) & 385 + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) ) & 386 + zbt1 * psxy(ji,jj) 387 psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 388 psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 389 END DO 390 END DO 391 392 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 393 DO ji = 1, jpi 394 zbt = zbet(ji,jj) 395 zbt1 = ( 1.0 - zbet(ji,jj) ) 396 psm(ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 397 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj) 398 zalf1 = 1.0 - zalf 399 ztemp = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 400 ps0(ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 401 psy(ji,jj) = zbt * psy(ji,jj) & 402 & + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 403 psyy(ji,jj) = zbt * psyy(ji,jj) & 404 & + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj) & 405 & + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 406 psxy(ji,jj) = zbt * psxy(ji,jj) & 407 & + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) & 408 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 409 psx(ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 410 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 411 END DO 412 END DO 325 ! Readjust moments remaining in the box. 326 DO jj = 2, jpj 327 DO ji = 1, jpi 328 zbt = zbet(ji,jj-1) 329 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 330 ! 331 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 332 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 333 psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 334 psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj) 335 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 336 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 337 psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 338 END DO 339 END DO 340 341 ! Put the temporary moments into appropriate neighboring boxes. 342 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 343 DO ji = 1, jpi 344 zbt = zbet(ji,jj-1) 345 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 346 psm(ji,jj) = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 347 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj) 348 zalf1 = 1.0 - zalf 349 ztemp = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 350 ! 351 ps0(ji,jj) = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj) 352 psy(ji,jj) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) & 353 & + zbt1 * psy(ji,jj) 354 psyy(ji,jj) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj) & 355 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 356 & + zbt1 * psyy(ji,jj) 357 psxy(ji,jj) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj) & 358 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) ) & 359 & + zbt1 * psxy(ji,jj) 360 psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 361 psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 362 END DO 363 END DO 364 ! 365 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 366 DO ji = 1, jpi 367 zbt = zbet(ji,jj) 368 zbt1 = ( 1.0 - zbet(ji,jj) ) 369 psm(ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 370 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj) 371 zalf1 = 1.0 - zalf 372 ztemp = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 373 ps0(ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 374 psy(ji,jj) = zbt * psy(ji,jj) & 375 & + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 376 psyy(ji,jj) = zbt * psyy(ji,jj) & 377 & + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj) & 378 & + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 379 psxy(ji,jj) = zbt * psxy(ji,jj) & 380 & + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) & 381 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 382 psx(ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 383 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 384 END DO 385 END DO 413 386 414 387 !-- Lateral boundary conditions 415 CALL lbc_lnk( psm , 'T', 1. ) 416 CALL lbc_lnk( ps0 , 'T', 1. ) 417 CALL lbc_lnk( psx , 'T', -1. ) 418 CALL lbc_lnk( psxx, 'T', 1. ) 419 CALL lbc_lnk( psy , 'T', -1. ) 420 CALL lbc_lnk( psyy, 'T', 1. ) 388 CALL lbc_lnk( psm , 'T', 1. ) ; CALL lbc_lnk( ps0 , 'T', 1. ) 389 CALL lbc_lnk( psx , 'T', -1. ) ; CALL lbc_lnk( psy , 'T', -1. ) ! caution gradient ==> the sign changes 390 CALL lbc_lnk( psxx, 'T', 1. ) ; CALL lbc_lnk( psyy, 'T', 1. ) 421 391 CALL lbc_lnk( psxy, 'T', 1. ) 422 392 … … 427 397 CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :') 428 398 ENDIF 429 399 ! 430 400 END SUBROUTINE lim_adv_y_2 431 401 … … 434 404 !! Default option Dummy module NO LIM 2.0 sea-ice model 435 405 !!---------------------------------------------------------------------- 436 CONTAINS437 SUBROUTINE lim_adv_x_2 ! Empty routine438 END SUBROUTINE lim_adv_x_2439 SUBROUTINE lim_adv_y_2 ! Empty routine440 END SUBROUTINE lim_adv_y_2441 442 406 #endif 443 407 408 !!====================================================================== 444 409 END MODULE limadv_2
Note: See TracChangeset
for help on using the changeset viewer.
