Changeset 1530 for trunk/NEMO/LIM_SRC_3/limadv.F90
- Timestamp:
- 2009-07-23T18:22:51+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limadv.F90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limadv.F90
r1529 r1530 4 4 !! LIM sea-ice model : sea-ice advection 5 5 !!====================================================================== 6 !! History : LIM ! 2008-03 (M. Vancoppenolle) LIM-3 from LIM-2 code 7 !! 3.2 ! 2009-06 (F. Dupont) correct a error in the North fold b. c. 8 !!-------------------------------------------------------------------- 6 9 #if defined key_lim3 7 10 !!---------------------------------------------------------------------- … … 11 14 !! lim_adv_y : advection of sea ice on y axis 12 15 !!---------------------------------------------------------------------- 13 !! * Modules used14 16 USE dom_oce 15 17 USE dom_ice 16 18 USE ice 19 USE lbclnk 17 20 USE in_out_manager ! I/O manager 18 USE lbclnk19 21 USE prtctl ! Print control 20 22 … … 22 24 PRIVATE 23 25 24 !! * Routine accessibility 25 PUBLIC lim_adv_x ! called by lim_trp 26 PUBLIC lim_adv_y ! called by lim_trp 26 PUBLIC lim_adv_x ! called by lim_trp 27 PUBLIC lim_adv_y ! called by lim_trp 28 29 REAL(wp) :: epsi20 = 1.e-20 ! constant values 30 REAL(wp) :: rzero = 0.e0 ! - - 31 REAL(wp) :: rone = 1.e0 ! - - 27 32 28 33 !! * Substitutions 29 34 # 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 3.0, UCL-LOCEAN-IPSL (2005) 35 !!---------------------------------------------------------------------- 36 !! NEMO/LIM 3.2, UCL-LOCEAN-IPSL (2009) 38 37 !! $Id$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 39 !!---------------------------------------------------------------------- 41 40 … … 48 47 !! 49 48 !! ** purpose : Computes and adds the advection trend to sea-ice 50 !! variable on xaxis49 !! variable on i-axis 51 50 !! 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 51 !! ** method : Uses Prather second order scheme that advects tracers 52 !! but also theirquadratic forms. The method preserves 53 !! tracer structures by conserving second order moments. 54 !! 55 !! Reference: Prather, 1986, JGR, 91, D6. 6671-6681. 63 56 !!-------------------------------------------------------------------- 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 ! " " 57 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for the time step 58 REAL(wp) , INTENT(in ) :: pcrh ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 59 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 60 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psm ! area 61 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ps0 ! field to be advected 62 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments 63 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 64 !! 65 INTEGER :: ji, jj ! dummy loop indices 66 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp, zin0 ! temporary scalars 67 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 68 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 69 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 70 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - 71 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 87 72 !--------------------------------------------------------------------- 88 73 89 74 ! Limitation of moments. 90 75 91 zrdt = rdt_ice * pdf! If ice drift field is too fast, use an appropriate time step for advection.76 zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 92 77 93 78 DO jj = 1, jpj … … 120 105 zalf1 = 1.0 - zalf 121 106 zalf1q = zalf1 * zalf1 122 zfm (ji,jj) = zalf * psm(ji,jj) 123 zf0 (ji,jj) = zalf * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) ) ) 124 zfx (ji,jj) = zalfq * ( psx(ji,jj) + 3.0 * zalf1 * psxx(ji,jj) ) 125 zfxx(ji,jj) = zalf * zalfq * psxx(ji,jj) 126 zfy (ji,jj) = zalf * ( psy(ji,jj) + zalf1 * psxy(ji,jj) ) 127 zfxy(ji,jj) = zalfq * psxy(ji,jj) 128 zfyy(ji,jj) = zalf * psyy(ji,jj) 107 ! 108 zfm (ji,jj) = zalf * psm (ji,jj) 109 zf0 (ji,jj) = zalf * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) ) ) 110 zfx (ji,jj) = zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) ) 111 zfxx(ji,jj) = zalf * psxx(ji,jj) * zalfq 112 zfy (ji,jj) = zalf * ( psy (ji,jj) + zalf1 * psxy(ji,jj) ) 113 zfxy(ji,jj) = zalfq * psxy(ji,jj) 114 zfyy(ji,jj) = zalf * psyy(ji,jj) 129 115 130 116 ! Readjust moments remaining in the box. … … 148 134 zalf1q = zalf1 * zalf1 149 135 zalg1q(ji,jj) = zalf1q 150 zfm (ji,jj) = zfm (ji,jj) + zalf * psm(ji+1,jj) 151 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0(ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) ) 152 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx(ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) ) 153 zfxx (ji,jj) = zfxx(ji,jj) + zalf * zalfq * psxx(ji+1,jj) 154 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy(ji+1,jj) - zalf1 * psxy(ji+1,jj) ) 155 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj) 156 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj) 136 ! 137 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj) 138 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) ) 139 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) ) 140 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj) * zalfq 141 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) ) 142 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj) 143 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj) 157 144 END DO 158 145 END DO … … 162 149 zbt = zbet(ji-1,jj) 163 150 zbt1 = 1.0 - zbet(ji-1,jj) 151 ! 164 152 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) ) 165 153 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) ) … … 181 169 zalf1 = 1.0 - zalf 182 170 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) & 171 ! 172 ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj) 173 psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 174 psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj) & 186 175 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 187 & + zbt1 * psxx(ji,jj)176 & + zbt1 * psxx(ji,jj) 188 177 psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj) & 189 178 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj) ) ) & 190 & + zbt1 * psxy(ji,jj)179 & + zbt1 * psxy(ji,jj) 191 180 psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj) 192 181 psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj) … … 201 190 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj) 202 191 zalf1 = 1.0 - zalf 203 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) ) ) 192 ztemp = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 193 ! 194 ps0(ji,jj) = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 195 psx(ji,jj) = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 196 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj) & 197 & + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) & 198 & + ( zalf1 - zalf ) * ztemp ) ) 199 psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) & 200 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) ) ) 213 201 psy(ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 214 202 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) … … 217 205 218 206 !-- 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. ) 207 CALL lbc_lnk( psm , 'T', 1. ) ; CALL lbc_lnk( ps0 , 'T', 1. ) 208 CALL lbc_lnk( psx , 'T', -1. ) ; CALL lbc_lnk( psy , 'T', -1. ) ! caution gradient ==> the sign changes 209 CALL lbc_lnk( psxx, 'T', 1. ) ; CALL lbc_lnk( psyy, 'T', 1. ) 225 210 CALL lbc_lnk( psxy, 'T', 1. ) 226 211 … … 231 216 CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :') 232 217 ENDIF 233 218 ! 234 219 END SUBROUTINE lim_adv_x 235 220 … … 241 226 !! 242 227 !! ** purpose : Computes and adds the advection trend to sea-ice 243 !! variable on y axis228 !! variable on y axis 244 229 !! 245 230 !! ** method : Uses Prather second order scheme that advects tracers 246 !! but also their quadratic forms. The method preserves tracer 247 !! structures by conserving second order moments. 231 !! but also their quadratic forms. The method preserves 232 !! tracer structures by conserving second order moments. 233 !! 234 !! Reference: Prather, 1986, JGR, 91, D6. 6671-6681. 235 !!--------------------------------------------------------------------- 236 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for the time step 237 REAL(wp) , INTENT(in ) :: pcrh ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 238 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 239 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psm ! area 240 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ps0 ! field to be advected 241 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments 242 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 248 243 !! 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 254 !!--------------------------------------------------------------------- 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 ! " " 244 INTEGER :: ji, jj ! dummy loop indices 245 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp, zin0 ! temporary scalars 246 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 247 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 248 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 249 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - 250 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 278 251 !--------------------------------------------------------------------- 279 252 … … 290 263 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) ) ) 291 264 zin0 = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask 265 ! 292 266 ps0 (ji,jj) = zslpmax 293 psx (ji,jj) 294 psxx(ji,jj) 267 psx (ji,jj) = psx (ji,jj) * zin0 268 psxx(ji,jj) = psxx(ji,jj) * zin0 295 269 psy (ji,jj) = zs1new * zin0 296 270 psyy(ji,jj) = zs2new * zin0 … … 300 274 301 275 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 302 psm (:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20)276 psm(:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 303 277 304 278 ! Calculate fluxes and moments between boxes j<-->j+1 305 DO jj = 1, jpj 306 DO ji = 1, jpi 307 ! Flux from j to j+1 WHEN v GT 0 279 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 280 DO ji = 1, jpi 308 281 zbet(ji,jj) = MAX( rzero, SIGN( rone, pvt(ji,jj) ) ) 309 282 zalf = MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) … … 311 284 zalf1 = 1.0 - zalf 312 285 zalf1q = zalf1 * zalf1 286 ! 313 287 zfm (ji,jj) = zalf * psm(ji,jj) 314 288 zf0 (ji,jj) = zalf * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj) + (zalf1-zalf) * psyy(ji,jj) ) ) … … 318 292 zfxy(ji,jj) = zalfq * psxy(ji,jj) 319 293 zfxx(ji,jj) = zalf * psxx(ji,jj) 320 294 ! 321 295 ! Readjust moments remaining in the box. 322 296 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj) … … 329 303 END DO 330 304 END DO 331 305 ! 332 306 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 333 307 DO ji = 1, jpi … … 339 313 zalf1q = zalf1 * zalf1 340 314 zalg1q(ji,jj) = zalf1q 341 zfm (ji,jj) = zfm (ji,jj) + zalf * psm(ji,jj+1) 342 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 343 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 344 zfyy (ji,jj) = zfyy(ji,jj) + zalf * zalfq * psyy(ji,jj+1) 345 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 346 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 347 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1) 315 ! 316 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1) 317 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 318 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 319 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1) * zalfq 320 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 321 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 322 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1) 348 323 END DO 349 324 END DO … … 354 329 zbt = zbet(ji,jj-1) 355 330 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 331 ! 356 332 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 357 333 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 358 334 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) 335 psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj) 360 336 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 361 337 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) … … 373 349 zalf1 = 1.0 - zalf 374 350 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 351 ! 352 ps0(ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj) 377 353 psy(ji,jj) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) & 378 & + zbt1 * psy(ji,jj) 379 354 & + zbt1 * psy(ji,jj) 380 355 psyy(ji,jj) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj) & 381 356 & + 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) 357 & + zbt1 * psyy(ji,jj) 358 psxy(ji,jj) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj) & 359 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) ) & 360 & + zbt1 * psxy(ji,jj) 387 361 psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 388 362 psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) … … 397 371 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj) 398 372 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) ) 373 ztemp = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj) 374 ps0 (ji,jj) = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 375 psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 376 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj) & 377 & + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) & 378 & + ( zalf1 - zalf ) * ztemp ) ) 379 psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) & 380 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 381 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 382 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 411 383 END DO 412 384 END DO 413 385 414 386 !-- 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. ) 387 CALL lbc_lnk( psm , 'T', 1. ) ; CALL lbc_lnk( ps0 , 'T', 1. ) 388 CALL lbc_lnk( psx , 'T', -1. ) ; CALL lbc_lnk( psy , 'T', -1. ) ! caution gradient ==> the sign changes 389 CALL lbc_lnk( psxx, 'T', 1. ) ; CALL lbc_lnk( psyy, 'T', 1. ) 421 390 CALL lbc_lnk( psxy, 'T', 1. ) 422 391 … … 427 396 CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :') 428 397 ENDIF 429 398 ! 430 399 END SUBROUTINE lim_adv_y 431 400 401 432 402 #else 433 403 !!---------------------------------------------------------------------- 434 404 !! Default option Dummy module NO LIM sea-ice model 435 405 !!---------------------------------------------------------------------- 436 CONTAINS437 SUBROUTINE lim_adv_x ! Empty routine438 END SUBROUTINE lim_adv_x439 SUBROUTINE lim_adv_y ! Empty routine440 END SUBROUTINE lim_adv_y441 442 406 #endif 443 407 408 !!====================================================================== 444 409 END MODULE limadv
Note: See TracChangeset
for help on using the changeset viewer.
