Changeset 6079
- Timestamp:
- 2015-12-17T11:08:49+01:00 (9 years ago)
- Location:
- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO
- Files:
-
- 5 deleted
- 56 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_2
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r5901 r6079 22 22 USE oce 23 23 USE dom_oce 24 USE sol_oce25 24 USE agrif_oce 26 25 USE phycst … … 29 28 USE lib_mpp 30 29 USE wrk_nemo 31 USE dynspg_oce32 30 USE zdf_oce 33 31 … … 38 36 39 37 PUBLIC Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts 40 PUBLIC interpun, interpvn , interpun2d, interpvn2d38 PUBLIC interpun, interpvn 41 39 PUBLIC interptsn, interpsshn 42 40 PUBLIC interpunb, interpvnb, interpub2b, interpvb2b … … 80 78 !! 81 79 INTEGER :: ji,jj,jk, j1,j2, i1,i2 82 REAL(wp) :: timeref 83 REAL(wp) :: z2dt, znugdt 84 REAL(wp) :: zrhox, zrhoy 85 REAL(wp), POINTER, DIMENSION(:,:) :: spgv1, spgu1 80 REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb 86 81 !!---------------------------------------------------------------------- 87 82 88 83 IF( Agrif_Root() ) RETURN 89 84 90 CALL wrk_alloc( jpi, jpj, spgv1, spgu1)85 CALL wrk_alloc( jpi, jpj, zub, zvb ) 91 86 92 87 Agrif_SpecialValue=0. … … 96 91 CALL Agrif_Bc_variable(vn_interp_id,procname=interpvn) 97 92 98 #if defined key_dynspg_flt99 CALL Agrif_Bc_variable(e1u_id,calledweight=1., procname=interpun2d)100 CALL Agrif_Bc_variable(e2v_id,calledweight=1., procname=interpvn2d)101 #endif102 103 93 Agrif_UseSpecialValue = .FALSE. 104 105 zrhox = Agrif_Rhox() 106 zrhoy = Agrif_Rhoy() 107 108 timeref = 1. 109 ! time step: leap-frog 110 z2dt = 2. * rdt 111 ! time step: Euler if restart from rest 112 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 113 ! coefficients 114 znugdt = grav * z2dt 115 94 116 95 ! prevent smoothing in ghost cells 117 96 i1=1 … … 126 105 127 106 IF((nbondi == -1).OR.(nbondi == 2)) THEN 128 #if defined key_dynspg_flt 129 DO jk=1,jpkm1 107 108 ! Smoothing 109 ! --------- 110 IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 111 ua_b(2,:)=0._wp 112 DO jk=1,jpkm1 113 DO jj=1,jpj 114 ua_b(2,jj) = ua_b(2,jj) + fse3u_a(2,jj,jk) * ua(2,jj,jk) 115 END DO 116 END DO 117 DO jj=1,jpj 118 ua_b(2,jj) = ua_b(2,jj) * hur_a(2,jj) 119 END DO 120 ENDIF 121 122 DO jk=1,jpkm1 ! Smooth 130 123 DO jj=j1,j2 131 ua(2,jj,jk) = (ua(2,jj,jk) - z2dt * znugdt * laplacu(2,jj))*umask(2,jj,jk)132 END DO133 END DO134 135 spgu(2,:)=0. 136 124 ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk)) 125 ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk) 126 END DO 127 END DO 128 129 zub(2,:)=0._wp ! Correct transport 137 130 DO jk=1,jpkm1 138 131 DO jj=1,jpj 139 spgu(2,jj)=spgu(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk) 140 END DO 141 END DO 142 132 zub(2,jj) = zub(2,jj) + fse3u_a(2,jj,jk) * ua(2,jj,jk) 133 END DO 134 END DO 143 135 DO jj=1,jpj 144 IF (umask(2,jj,1).NE.0.) THEN 145 spgu(2,jj)=spgu(2,jj)/hu(2,jj) 146 ENDIF 147 END DO 148 #else 149 spgu(2,:) = ua_b(2,:) 150 #endif 151 152 DO jk=1,jpkm1 136 zub(2,jj) = zub(2,jj) * hur_a(2,jj) 137 END DO 138 139 DO jk=1,jpkm1 140 DO jj=1,jpj 141 ua(2,jj,jk) = (ua(2,jj,jk)+ua_b(2,jj)-zub(2,jj))*umask(2,jj,jk) 142 END DO 143 END DO 144 145 ! Set tangential velocities to time splitting estimate 146 !----------------------------------------------------- 147 IF ( ln_dynspg_ts) THEN 148 zvb(2,:)=0._wp 149 DO jk=1,jpkm1 150 DO jj=1,jpj 151 zvb(2,jj) = zvb(2,jj) + fse3v_a(2,jj,jk) * va(2,jj,jk) 152 END DO 153 END DO 154 DO jj=1,jpj 155 zvb(2,jj) = zvb(2,jj) * hvr_a(2,jj) 156 END DO 157 DO jk=1,jpkm1 158 DO jj=1,jpj 159 va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-zvb(2,jj))*vmask(2,jj,jk) 160 END DO 161 END DO 162 ENDIF 163 164 ! Mask domain edges: 165 !------------------- 166 DO jk=1,jpkm1 167 DO jj=1,jpj 168 ua(1,jj,jk) = 0._wp 169 va(1,jj,jk) = 0._wp 170 END DO 171 END DO 172 173 ENDIF 174 175 IF((nbondi == 1).OR.(nbondi == 2)) THEN 176 177 ! Smoothing 178 ! --------- 179 IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 180 ua_b(nlci-2,:)=0._wp 181 DO jk=1,jpkm1 182 DO jj=1,jpj 183 ua_b(nlci-2,jj) = ua_b(nlci-2,jj) + fse3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk) 184 END DO 185 END DO 186 DO jj=1,jpj 187 ua_b(nlci-2,jj) = ua_b(nlci-2,jj) * hur_a(nlci-2,jj) 188 END DO 189 ENDIF 190 191 DO jk=1,jpkm1 ! Smooth 153 192 DO jj=j1,j2 154 ua(2,jj,jk) = 0.25*(ua(1,jj,jk)+2.*ua(2,jj,jk)+ua(3,jj,jk)) 155 ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk) 156 END DO 157 END DO 158 159 spgu1(2,:)=0. 160 193 ua(nlci-2,jj,jk) = 0.25_wp*(ua(nlci-3,jj,jk)+2._wp*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk)) 194 ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk) 195 END DO 196 END DO 197 198 zub(nlci-2,:)=0._wp ! Correct transport 161 199 DO jk=1,jpkm1 162 200 DO jj=1,jpj 163 spgu1(2,jj)=spgu1(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk) 164 END DO 165 END DO 166 201 zub(nlci-2,jj) = zub(nlci-2,jj) + fse3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk) 202 END DO 203 END DO 167 204 DO jj=1,jpj 168 IF (umask(2,jj,1).NE.0.) THEN 169 spgu1(2,jj)=spgu1(2,jj)/hu(2,jj) 170 ENDIF 171 END DO 172 173 DO jk=1,jpkm1 174 DO jj=j1,j2 175 ua(2,jj,jk) = (ua(2,jj,jk)+spgu(2,jj)-spgu1(2,jj))*umask(2,jj,jk) 176 END DO 177 END DO 178 179 #if defined key_dynspg_ts 205 zub(nlci-2,jj) = zub(nlci-2,jj) * hur_a(nlci-2,jj) 206 END DO 207 208 DO jk=1,jpkm1 209 DO jj=1,jpj 210 ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)+ua_b(nlci-2,jj)-zub(nlci-2,jj))*umask(nlci-2,jj,jk) 211 END DO 212 END DO 213 180 214 ! Set tangential velocities to time splitting estimate 181 spgv1(2,:)=0. 182 DO jk=1,jpkm1 215 !----------------------------------------------------- 216 IF ( ln_dynspg_ts) THEN 217 zvb(nlci-1,:)=0._wp 218 DO jk=1,jpkm1 219 DO jj=1,jpj 220 zvb(nlci-1,jj) = zvb(nlci-1,jj) + fse3v_a(nlci-1,jj,jk) * va(nlci-1,jj,jk) 221 END DO 222 END DO 183 223 DO jj=1,jpj 184 spgv1(2,jj)=spgv1(2,jj)+fse3v_a(2,jj,jk)*va(2,jj,jk) 185 END DO 186 END DO 187 DO jj=1,jpj 188 spgv1(2,jj)=spgv1(2,jj)*hvr_a(2,jj) 189 END DO 224 zvb(nlci-1,jj) = zvb(nlci-1,jj) * hvr_a(nlci-1,jj) 225 END DO 226 DO jk=1,jpkm1 227 DO jj=1,jpj 228 va(nlci-1,jj,jk) = (va(nlci-1,jj,jk)+va_b(nlci-1,jj)-zvb(nlci-1,jj))*vmask(nlci-1,jj,jk) 229 END DO 230 END DO 231 ENDIF 232 233 ! Mask domain edges: 234 !------------------- 190 235 DO jk=1,jpkm1 191 236 DO jj=1,jpj 192 va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-spgv1(2,jj))*vmask(2,jj,jk) 193 END DO 194 END DO 195 #endif 196 197 ENDIF 198 199 IF((nbondi == 1).OR.(nbondi == 2)) THEN 200 #if defined key_dynspg_flt 201 DO jk=1,jpkm1 202 DO jj=j1,j2 203 ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)- z2dt * znugdt * laplacu(nlci-2,jj))*umask(nlci-2,jj,jk) 204 END DO 205 END DO 206 spgu(nlci-2,:)=0. 207 DO jk=1,jpkm1 208 DO jj=1,jpj 209 spgu(nlci-2,jj)=spgu(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk) 210 ENDDO 211 ENDDO 212 DO jj=1,jpj 213 IF (umask(nlci-2,jj,1).NE.0.) THEN 214 spgu(nlci-2,jj)=spgu(nlci-2,jj)/hu(nlci-2,jj) 215 ENDIF 216 END DO 217 #else 218 spgu(nlci-2,:) = ua_b(nlci-2,:) 219 #endif 220 DO jk=1,jpkm1 221 DO jj=j1,j2 222 ua(nlci-2,jj,jk) = 0.25*(ua(nlci-3,jj,jk)+2.*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk)) 223 224 ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk) 225 226 END DO 227 END DO 228 spgu1(nlci-2,:)=0. 229 DO jk=1,jpkm1 230 DO jj=1,jpj 231 spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)*umask(nlci-2,jj,jk) 232 END DO 233 END DO 234 DO jj=1,jpj 235 IF (umask(nlci-2,jj,1).NE.0.) THEN 236 spgu1(nlci-2,jj)=spgu1(nlci-2,jj)/hu(nlci-2,jj) 237 ENDIF 238 END DO 239 DO jk=1,jpkm1 240 DO jj=j1,j2 241 ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)+spgu(nlci-2,jj)-spgu1(nlci-2,jj))*umask(nlci-2,jj,jk) 242 END DO 243 END DO 244 245 #if defined key_dynspg_ts 237 ua(nlci-1,jj,jk) = 0._wp 238 va(nlci ,jj,jk) = 0._wp 239 END DO 240 END DO 241 242 ENDIF 243 244 IF((nbondj == -1).OR.(nbondj == 2)) THEN 245 246 ! Smoothing 247 ! --------- 248 IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 249 va_b(:,2)=0._wp 250 DO jk=1,jpkm1 251 DO ji=1,jpi 252 va_b(ji,2) = va_b(ji,2) + fse3v_a(ji,2,jk) * va(ji,2,jk) 253 END DO 254 END DO 255 DO ji=1,jpi 256 va_b(ji,2) = va_b(ji,2) * hvr_a(ji,2) 257 END DO 258 ENDIF 259 260 DO jk=1,jpkm1 ! Smooth 261 DO ji=i1,i2 262 va(ji,2,jk)=0.25_wp*(va(ji,1,jk)+2._wp*va(ji,2,jk)+va(ji,3,jk)) 263 va(ji,2,jk)=va(ji,2,jk)*vmask(ji,2,jk) 264 END DO 265 END DO 266 267 zvb(:,2)=0._wp ! Correct transport 268 DO jk=1,jpkm1 269 DO ji=1,jpi 270 zvb(ji,2) = zvb(ji,2) + fse3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk) 271 END DO 272 END DO 273 DO ji=1,jpi 274 zvb(ji,2) = zvb(ji,2) * hvr_a(ji,2) 275 END DO 276 DO jk=1,jpkm1 277 DO ji=1,jpi 278 va(ji,2,jk) = (va(ji,2,jk)+va_b(ji,2)-zvb(ji,2))*vmask(ji,2,jk) 279 END DO 280 END DO 281 246 282 ! Set tangential velocities to time splitting estimate 247 spgv1(nlci-1,:)=0._wp 248 DO jk=1,jpkm1 249 DO jj=1,jpj 250 spgv1(nlci-1,jj)=spgv1(nlci-1,jj)+fse3v_a(nlci-1,jj,jk)*va(nlci-1,jj,jk)*vmask(nlci-1,jj,jk) 251 END DO 252 END DO 253 254 DO jj=1,jpj 255 spgv1(nlci-1,jj)=spgv1(nlci-1,jj)*hvr_a(nlci-1,jj) 256 END DO 257 258 DO jk=1,jpkm1 259 DO jj=1,jpj 260 va(nlci-1,jj,jk) = (va(nlci-1,jj,jk)+va_b(nlci-1,jj)-spgv1(nlci-1,jj))*vmask(nlci-1,jj,jk) 261 END DO 262 END DO 263 #endif 264 265 ENDIF 266 267 IF((nbondj == -1).OR.(nbondj == 2)) THEN 268 269 #if defined key_dynspg_flt 270 DO jk=1,jpkm1 283 !----------------------------------------------------- 284 IF ( ln_dynspg_ts ) THEN 285 zub(:,2)=0._wp 286 DO jk=1,jpkm1 287 DO ji=1,jpi 288 zub(ji,2) = zub(ji,2) + fse3u_a(ji,2,jk) * ua(ji,2,jk) * umask(ji,2,jk) 289 END DO 290 END DO 271 291 DO ji=1,jpi 272 va(ji,2,jk) = (va(ji,2,jk) - z2dt * znugdt * laplacv(ji,2))*vmask(ji,2,jk) 273 END DO 274 END DO 275 276 spgv(:,2)=0. 277 292 zub(ji,2) = zub(ji,2) * hur_a(ji,2) 293 END DO 294 295 DO jk=1,jpkm1 296 DO ji=1,jpi 297 ua(ji,2,jk) = (ua(ji,2,jk)+ua_b(ji,2)-zub(ji,2))*umask(ji,2,jk) 298 END DO 299 END DO 300 ENDIF 301 302 ! Mask domain edges: 303 !------------------- 278 304 DO jk=1,jpkm1 279 305 DO ji=1,jpi 280 spgv(ji,2)=spgv(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk) 281 END DO 282 END DO 283 306 ua(ji,1,jk) = 0._wp 307 va(ji,1,jk) = 0._wp 308 END DO 309 END DO 310 311 ENDIF 312 313 IF((nbondj == 1).OR.(nbondj == 2)) THEN 314 ! Smoothing 315 ! --------- 316 IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 317 va_b(:,nlcj-2)=0._wp 318 DO jk=1,jpkm1 319 DO ji=1,jpi 320 va_b(ji,nlcj-2) = va_b(ji,nlcj-2) + fse3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) 321 END DO 322 END DO 323 DO ji=1,jpi 324 va_b(ji,nlcj-2) = va_b(ji,nlcj-2) * hvr_a(ji,nlcj-2) 325 END DO 326 ENDIF 327 328 DO jk=1,jpkm1 ! Smooth 329 DO ji=i1,i2 330 va(ji,nlcj-2,jk)=0.25_wp*(va(ji,nlcj-3,jk)+2._wp*va(ji,nlcj-2,jk)+va(ji,nlcj-1,jk)) 331 va(ji,nlcj-2,jk)=va(ji,nlcj-2,jk)*vmask(ji,nlcj-2,jk) 332 END DO 333 END DO 334 335 zvb(:,nlcj-2)=0._wp ! Correct transport 336 DO jk=1,jpkm1 337 DO ji=1,jpi 338 zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + fse3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk) 339 END DO 340 END DO 284 341 DO ji=1,jpi 285 IF (vmask(ji,2,1).NE.0.) THEN 286 spgv(ji,2)=spgv(ji,2)/hv(ji,2) 287 ENDIF 288 END DO 289 #else 290 spgv(:,2)=va_b(:,2) 291 #endif 292 293 DO jk=1,jpkm1 294 DO ji=i1,i2 295 va(ji,2,jk)=0.25*(va(ji,1,jk)+2.*va(ji,2,jk)+va(ji,3,jk)) 296 va(ji,2,jk)=va(ji,2,jk)*vmask(ji,2,jk) 297 END DO 298 END DO 299 300 spgv1(:,2)=0. 301 342 zvb(ji,nlcj-2) = zvb(ji,nlcj-2) * hvr_a(ji,nlcj-2) 343 END DO 302 344 DO jk=1,jpkm1 303 345 DO ji=1,jpi 304 spgv1(ji,2)=spgv1(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)*vmask(ji,2,jk) 305 END DO 306 END DO 307 308 DO ji=1,jpi 309 IF (vmask(ji,2,1).NE.0.) THEN 310 spgv1(ji,2)=spgv1(ji,2)/hv(ji,2) 311 ENDIF 312 END DO 313 314 DO jk=1,jpkm1 346 va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)+va_b(ji,nlcj-2)-zvb(ji,nlcj-2))*vmask(ji,nlcj-2,jk) 347 END DO 348 END DO 349 350 ! Set tangential velocities to time splitting estimate 351 !----------------------------------------------------- 352 IF ( ln_dynspg_ts ) THEN 353 zub(:,nlcj-1)=0._wp 354 DO jk=1,jpkm1 355 DO ji=1,jpi 356 zub(ji,nlcj-1) = zub(ji,nlcj-1) + fse3u_a(ji,nlcj-1,jk) * ua(ji,nlcj-1,jk) * umask(ji,nlcj-1,jk) 357 END DO 358 END DO 315 359 DO ji=1,jpi 316 va(ji,2,jk) = (va(ji,2,jk)+spgv(ji,2)-spgv1(ji,2))*vmask(ji,2,jk) 317 END DO 318 END DO 319 320 #if defined key_dynspg_ts 321 ! Set tangential velocities to time splitting estimate 322 spgu1(:,2)=0._wp 360 zub(ji,nlcj-1) = zub(ji,nlcj-1) * hur_a(ji,nlcj-1) 361 END DO 362 363 DO jk=1,jpkm1 364 DO ji=1,jpi 365 ua(ji,nlcj-1,jk) = (ua(ji,nlcj-1,jk)+ua_b(ji,nlcj-1)-zub(ji,nlcj-1))*umask(ji,nlcj-1,jk) 366 END DO 367 END DO 368 ENDIF 369 370 ! Mask domain edges: 371 !------------------- 323 372 DO jk=1,jpkm1 324 373 DO ji=1,jpi 325 spgu1(ji,2)=spgu1(ji,2)+fse3u_a(ji,2,jk)*ua(ji,2,jk)*umask(ji,2,jk) 326 END DO 327 END DO 328 329 DO ji=1,jpi 330 spgu1(ji,2)=spgu1(ji,2)*hur_a(ji,2) 331 END DO 332 333 DO jk=1,jpkm1 334 DO ji=1,jpi 335 ua(ji,2,jk) = (ua(ji,2,jk)+ua_b(ji,2)-spgu1(ji,2))*umask(ji,2,jk) 336 END DO 337 END DO 338 #endif 339 ENDIF 340 341 IF((nbondj == 1).OR.(nbondj == 2)) THEN 342 343 #if defined key_dynspg_flt 344 DO jk=1,jpkm1 345 DO ji=1,jpi 346 va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)-z2dt * znugdt * laplacv(ji,nlcj-2))*vmask(ji,nlcj-2,jk) 347 END DO 348 END DO 349 350 351 spgv(:,nlcj-2)=0. 352 353 DO jk=1,jpkm1 354 DO ji=1,jpi 355 spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk) 356 END DO 357 END DO 358 359 DO ji=1,jpi 360 IF (vmask(ji,nlcj-2,1).NE.0.) THEN 361 spgv(ji,nlcj-2)=spgv(ji,nlcj-2)/hv(ji,nlcj-2) 362 ENDIF 363 END DO 364 365 #else 366 spgv(:,nlcj-2)=va_b(:,nlcj-2) 367 #endif 368 369 DO jk=1,jpkm1 370 DO ji=i1,i2 371 va(ji,nlcj-2,jk)=0.25*(va(ji,nlcj-3,jk)+2.*va(ji,nlcj-2,jk)+va(ji,nlcj-1,jk)) 372 va(ji,nlcj-2,jk) = va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk) 373 END DO 374 END DO 375 376 spgv1(:,nlcj-2)=0. 377 378 DO jk=1,jpkm1 379 DO ji=1,jpi 380 spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk) 381 END DO 382 END DO 383 384 DO ji=1,jpi 385 IF (vmask(ji,nlcj-2,1).NE.0.) THEN 386 spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)/hv(ji,nlcj-2) 387 ENDIF 388 END DO 389 390 DO jk=1,jpkm1 391 DO ji=1,jpi 392 va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)+spgv(ji,nlcj-2)-spgv1(ji,nlcj-2))*vmask(ji,nlcj-2,jk) 393 END DO 394 END DO 395 396 #if defined key_dynspg_ts 397 ! Set tangential velocities to time splitting estimate 398 spgu1(:,nlcj-1)=0._wp 399 DO jk=1,jpkm1 400 DO ji=1,jpi 401 spgu1(ji,nlcj-1)=spgu1(ji,nlcj-1)+fse3u_a(ji,nlcj-1,jk)*ua(ji,nlcj-1,jk) 402 END DO 403 END DO 404 405 DO ji=1,jpi 406 spgu1(ji,nlcj-1)=spgu1(ji,nlcj-1)*hur_a(ji,nlcj-1) 407 END DO 408 409 DO jk=1,jpkm1 410 DO ji=1,jpi 411 ua(ji,nlcj-1,jk) = (ua(ji,nlcj-1,jk)+ua_b(ji,nlcj-1)-spgu1(ji,nlcj-1))*umask(ji,nlcj-1,jk) 412 END DO 413 END DO 414 #endif 415 416 ENDIF 417 ! 418 CALL wrk_dealloc( jpi, jpj, spgv1, spgu1 ) 374 ua(ji,nlcj ,jk) = 0._wp 375 va(ji,nlcj-1,jk) = 0._wp 376 END DO 377 END DO 378 379 ENDIF 380 ! 381 CALL wrk_dealloc( jpi, jpj, zub, zvb ) 419 382 ! 420 383 END SUBROUTINE Agrif_dyn … … 687 650 END DO 688 651 END DO 652 tsa(nlci,j1:j2,k1:k2,jn) = 0._wp 689 653 ENDDO 690 654 ENDIF … … 706 670 END DO 707 671 END DO 672 tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp 708 673 ENDDO 709 674 ENDIF … … 724 689 END DO 725 690 END DO 691 tsa(1,j1:j2,k1:k2,jn) = 0._wp 726 692 END DO 727 693 ENDIF … … 742 708 END DO 743 709 END DO 710 tsa(i1:i2,1,k1:k2,jn) = 0._wp 744 711 ENDDO 745 712 ENDIF … … 828 795 END SUBROUTINE interpun 829 796 830 831 SUBROUTINE interpun2d(ptab,i1,i2,j1,j2,before)832 !!---------------------------------------------833 !! *** ROUTINE interpun ***834 !!---------------------------------------------835 !836 INTEGER, INTENT(in) :: i1,i2,j1,j2837 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab838 LOGICAL, INTENT(in) :: before839 !840 INTEGER :: ji,jj841 REAL(wp) :: ztref842 REAL(wp) :: zrhoy843 !!---------------------------------------------844 !845 ztref = 1.846 847 IF (before) THEN848 DO jj=j1,j2849 DO ji=i1,MIN(i2,nlci-1)850 ptab(ji,jj) = e2u(ji,jj) * ((gcx(ji+1,jj) - gcx(ji,jj))/e1u(ji,jj))851 END DO852 END DO853 ELSE854 zrhoy = Agrif_Rhoy()855 DO jj=j1,j2856 laplacu(i1:i2,jj) = ztref * (ptab(i1:i2,jj)/(zrhoy*e2u(i1:i2,jj))) !*umask(i1:i2,jj,1)857 END DO858 ENDIF859 !860 END SUBROUTINE interpun2d861 862 863 797 SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before) 864 798 !!--------------------------------------------- … … 895 829 ! 896 830 END SUBROUTINE interpvn 897 898 SUBROUTINE interpvn2d(ptab,i1,i2,j1,j2,before)899 !!---------------------------------------------900 !! *** ROUTINE interpvn ***901 !!---------------------------------------------902 !903 INTEGER, INTENT(in) :: i1,i2,j1,j2904 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab905 LOGICAL, INTENT(in) :: before906 !907 INTEGER :: ji,jj908 REAL(wp) :: zrhox909 REAL(wp) :: ztref910 !!---------------------------------------------911 !912 ztref = 1.913 IF (before) THEN914 !interpv entre 1 et k2 et interpv2d en jpkp1915 DO jj=j1,MIN(j2,nlcj-1)916 DO ji=i1,i2917 ptab(ji,jj) = e1v(ji,jj) * ((gcx(ji,jj+1) - gcx(ji,jj))/e2v(ji,jj)) * vmask(ji,jj,1)918 END DO919 END DO920 ELSE921 zrhox = Agrif_Rhox()922 DO ji=i1,i2923 laplacv(ji,j1:j2) = ztref * (ptab(ji,j1:j2)/(zrhox*e1v(ji,j1:j2)))924 END DO925 ENDIF926 !927 END SUBROUTINE interpvn2d928 831 929 832 SUBROUTINE interpunb(ptab,i1,i2,j1,j2,before,nb,ndir) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
r5901 r6079 11 11 USE lib_mpp 12 12 USE wrk_nemo 13 USE dynspg_oce14 13 USE zdf_oce ! vertical physics: ocean variables 15 14 … … 107 106 # endif 108 107 109 # if defined key_dynspg_ts 110 IF (ln_bt_fw) THEN 108 IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 111 109 ! Update time integrated transports 112 110 IF (mod(nbcline,nbclineupdate) == 0) THEN … … 128 126 ENDIF 129 127 END IF 130 # endif131 128 ! 132 129 nbcline = nbcline + 1 … … 360 357 ! 361 358 ! Update barotropic velocities: 362 #if defined key_dynspg_ts 363 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part364 zcorr = tabres(ji,jj) - un_b(ji,jj)365 ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)366 END IF367 #endif359 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 360 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 361 zcorr = tabres(ji,jj) - un_b(ji,jj) 362 ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 363 END IF 364 ENDIF 368 365 un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1) 369 366 ! … … 425 422 ! 426 423 ! Update barotropic velocities: 427 #if defined key_dynspg_ts 428 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part429 zcorr = tabres(ji,jj) - vn_b(ji,jj)430 vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)431 END IF432 #endif424 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 425 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 426 zcorr = tabres(ji,jj) - vn_b(ji,jj) 427 vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 428 END IF 429 ENDIF 433 430 vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1) 434 431 ! … … 470 467 END DO 471 468 ELSE 472 #if ! defined key_dynspg_ts 473 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 474 DO jj=j1,j2 475 DO ji=i1,i2 476 sshb(ji,jj) = sshb(ji,jj) & 477 & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 478 END DO 479 END DO 469 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 470 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 471 DO jj=j1,j2 472 DO ji=i1,i2 473 sshb(ji,jj) = sshb(ji,jj) & 474 & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 475 END DO 476 END DO 477 ENDIF 480 478 ENDIF 481 #endif 479 482 480 DO jj=j1,j2 483 481 DO ji=i1,i2 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90
r5901 r6079 4 4 USE oce 5 5 USE dom_oce 6 USE sol_oce7 6 USE agrif_oce 8 7 USE agrif_top_sponge -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r5901 r6079 104 104 USE dom_oce 105 105 USE nemogcm 106 USE sol_oce107 106 USE in_out_manager 108 107 USE agrif_opa_update … … 172 171 USE dom_oce 173 172 USE nemogcm 174 USE sol_oce175 173 USE lib_mpp 176 174 USE in_out_manager … … 210 208 CALL Agrif_Bc_variable(vn_sponge_id,calledweight=1.,procname=interpvn_sponge) 211 209 212 #if defined key_dynspg_ts213 210 Agrif_UseSpecialValue = .TRUE. 214 211 CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn ) 215 212 216 Agrif_UseSpecialValue = ln_spc_dyn 217 CALL Agrif_Bc_variable(unb_id,calledweight=1.,procname=interpunb) 218 CALL Agrif_Bc_variable(vnb_id,calledweight=1.,procname=interpvnb) 219 CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 220 CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 221 ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 ; hbdy_w(:) =0.e0 222 ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 ; hbdy_e(:) =0.e0 223 ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 ; hbdy_n(:) =0.e0 224 ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 ; hbdy_s(:) =0.e0 225 #endif 213 IF ( ln_dynspg_ts ) THEN 214 Agrif_UseSpecialValue = ln_spc_dyn 215 CALL Agrif_Bc_variable(unb_id,calledweight=1.,procname=interpunb) 216 CALL Agrif_Bc_variable(vnb_id,calledweight=1.,procname=interpvnb) 217 CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 218 CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 219 ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 ; hbdy_w(:) =0.e0 220 ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 ; hbdy_e(:) =0.e0 221 ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 ; hbdy_n(:) =0.e0 222 ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 ; hbdy_s(:) =0.e0 223 ENDIF 226 224 227 225 Agrif_UseSpecialValue = .FALSE. … … 278 276 ENDIF 279 277 ENDIF 278 279 ! Check free surface scheme 280 IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.& 281 & ( Agrif_Parent(ln_dynspg_exp).AND.ln_dynspg_ts ) ) THEN 282 WRITE(*,*) 'incompatible free surface scheme between grids' 283 WRITE(*,*) 'parent grid ln_dynspg_ts :', Agrif_Parent(ln_dynspg_ts ) 284 WRITE(*,*) 'parent grid ln_dynspg_exp :', Agrif_Parent(ln_dynspg_exp) 285 WRITE(*,*) 'child grid ln_dynspg_ts :', ln_dynspg_ts 286 WRITE(*,*) 'child grid ln_dynspg_exp :', ln_dynspg_exp 287 WRITE(*,*) 'those logicals should be identical' 288 STOP 289 ENDIF 290 280 291 ! check if masks and bathymetries match 281 292 IF(ln_chk_bathy) THEN … … 319 330 ! 320 331 Agrif_UseSpecialValueInUpdate = .FALSE. 321 nbcline = 0 332 nbcline = 0 322 333 lk_agrif_doupd = .FALSE. 323 334 ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OFF_SRC
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
r5901 r6079 107 107 USE oce, ONLY: zts => tsa 108 108 USE oce, ONLY: zuslp => ua , zvslp => va 109 USE oce, ONLY: zwslpi => ua_sv , zwslpj => va_sv109 USE zdf_oce, ONLY: zwslpi => avmu , zwslpj => avmv 110 110 USE oce, ONLY: zu => ub , zv => vb, zw => rke 111 111 ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OOO_SRC
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
r5901 r6079 34 34 USE zdfmxl ! Mixed layer depth 35 35 USE dom_oce, ONLY : ndastp 36 USE sol_oce, ONLY : gcx ! Solver variables defined in memory37 36 USE in_out_manager ! I/O manager 38 37 USE iom ! I/O module … … 114 113 CALL iom_rstput( kt, nitbkg_r, inum, 'en' , en ) 115 114 #endif 116 CALL iom_rstput( kt, nitbkg_r, inum, 'gcx' , gcx )117 115 ! 118 116 CALL iom_close( inum ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r5901 r6079 86 86 ! 87 87 INTEGER :: nb_bdy !: number of open boundary sets 88 INTEGER ,:: nb_jpk_bdy !: number of levels in the bdy data (set < 0 if consistent with planned run)88 INTEGER :: nb_jpk_bdy !: number of levels in the bdy data (set < 0 if consistent with planned run) 89 89 INTEGER, DIMENSION(jp_bdy) :: nn_rimwidth !: boundary rim width for Flow Relaxation Scheme 90 90 INTEGER :: nn_volctl !: = 0 the total volume will have the variability of the surface Flux E-P … … 131 131 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global !: workspace for reading in global data arrays (unstr. bdy) 132 132 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global_z !: workspace for reading in global depth arrays (unstr. bdy) 133 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global_dz 133 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global_dz !: workspace for reading in global depth arrays (unstr. bdy) 134 134 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2 !: workspace for reading in global data arrays (struct. bdy) 135 135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2_z !: workspace for reading in global depth arrays (struct. bdy) 136 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2_dz 136 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2_dz !: workspace for reading in global depth arrays (struct. bdy) 137 137 !$AGRIF_DO_NOT_TREAT 138 138 TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET :: idx_bdy !: bdy indices (local process) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r5620 r6079 29 29 USE iom ! IOM library 30 30 USE in_out_manager ! I/O logical units 31 USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag32 31 #if defined key_lim2 33 32 USE ice_2 … … 393 392 END DO ! ib_bdy 394 393 395 ! bg jchanut tschanges396 394 #if defined key_tide 397 ! Add tides if not split-explicit free surface else this is done in ts loop 398 IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 399 #endif 400 ! end jchanut tschanges 395 IF (ln_dynspg_ts) THEN ! Fill temporary arrays with slow-varying bdy data 396 DO ib_bdy = 1, nb_bdy ! Tidal component added in ts loop 397 IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 398 nblen => idx_bdy(ib_bdy)%nblen 399 nblenrim => idx_bdy(ib_bdy)%nblenrim 400 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 401 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen1(1)) = dta_bdy(ib_bdy)%ssh(1:ilen1(1)) 402 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen1(2)) = dta_bdy(ib_bdy)%u2d(1:ilen1(2)) 403 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen1(3)) = dta_bdy(ib_bdy)%v2d(1:ilen1(3)) 404 ENDIF 405 END DO 406 ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 407 ! 408 CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 409 ENDIF 410 #endif 401 411 402 412 IF ( ln_apr_obc ) THEN … … 428 438 !! 429 439 !!---------------------------------------------------------------------- 430 USE dynspg_oce, ONLY: lk_dynspg_ts431 440 !! 432 441 INTEGER :: ib_bdy, jfld, jstart, jend, ierror ! local indices … … 435 444 CHARACTER(len=100) :: cn_dir ! Root directory for location of data files 436 445 CHARACTER(len=100), DIMENSION(nb_bdy) :: cn_dir_array ! Root directory for location of data files 446 CHARACTER(len = 256):: clname ! temporary file name 437 447 LOGICAL :: ln_full_vel ! =T => full velocities in 3D boundary data 438 448 ! =F => baroclinic velocities in 3D boundary data … … 675 685 ! sea ice 676 686 IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 677 678 687 ! Test for types of ice input (lim2 or lim3) 679 CALL iom_open ( bn_a_i%clname, inum ) 680 id1 = iom_varid ( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 688 ! Build file name to find dimensions 689 clname=TRIM(bn_a_i%clname) 690 IF( .NOT. bn_a_i%ln_clim ) THEN 691 WRITE(clname, '(a,"_y",i4.4)' ) TRIM( bn_a_i%clname ), nyear ! add year 692 IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), nmonth ! add month 693 ELSE 694 IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( bn_a_i%clname ), nmonth ! add month 695 ENDIF 696 IF( bn_a_i%cltype == 'daily' .OR. bn_a_i%cltype(1:4) == 'week' ) & 697 & WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), nday ! add day 698 ! 699 CALL iom_open ( clname, inum ) 700 id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 681 701 CALL iom_close ( inum ) 682 !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. ) 683 !CALL iom_open ( bn_a_i%clname, inum ) 684 !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 702 685 703 IF ( zndims == 4 ) THEN 686 704 ll_bdylim3 = .TRUE. ! lim3 input -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
r4792 r6079 24 24 USE oce ! ocean dynamics and tracers 25 25 USE dom_oce ! ocean space and time domain 26 USE dynspg_oce27 26 USE bdy_oce ! ocean open boundary conditions 28 27 USE bdydyn2d ! open boundary conditions for barotropic solution … … 35 34 PRIVATE 36 35 37 PUBLIC bdy_dyn ! routine called in dynspg_flt (if lk_dynspg_flt) or 38 ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp) 36 PUBLIC bdy_dyn ! routine called in dyn_nxt 39 37 40 38 # include "domzgr_substitute.h90" -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
r5620 r6079 23 23 USE bdy_oce ! ocean open boundary conditions 24 24 USE bdylib ! BDY library routines 25 USE dynspg_oce ! for barotropic variables26 25 USE phycst ! physical constants 27 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r5620 r6079 33 33 ! USE tide_mod ! Useless ?? 34 34 USE fldread 35 USE dynspg_oce, ONLY: lk_dynspg_ts36 35 37 36 IMPLICIT NONE … … 54 53 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data 55 54 !$AGRIF_END_DO_NOT_TREAT 56 TYPE(OBC_DATA) , P RIVATE, DIMENSION(jp_bdy) :: dta_bdy_s !: bdy external data (slow component)55 TYPE(OBC_DATA) , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s !: bdy external data (slow component) 57 56 58 57 !!---------------------------------------------------------------------- … … 270 269 ENDIF 271 270 ! 272 IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during 273 ! time splitting integration 274 ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 275 ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 276 ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 277 dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 278 dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 279 dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 280 ENDIF 271 ! Allocate slow varying data in the case of time splitting: 272 ! Do it anyway because at this stage knowledge of free surface scheme is unknown 273 ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 274 ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 275 ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 276 dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 277 dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 278 dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 281 279 ! 282 280 ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 … … 397 395 !! 398 396 LOGICAL :: lk_first_btstp ! =.TRUE. if time splitting and first barotropic step 399 INTEGER, 397 INTEGER, DIMENSION(jpbgrd) :: ilen0 400 398 INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim ! short cuts 401 399 INTEGER :: itide, ib_bdy, ib, igrd ! loop indices … … 416 414 ! Absolute time from model initialization: 417 415 IF( PRESENT(kit) ) THEN 418 z_arg = ( kt + (kit+ 0.5_wp*(time_add-1)) / REAL(nn_baro,wp) ) * rdt416 z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt 419 417 ELSE 420 418 z_arg = ( kt + time_add ) * rdt … … 456 454 zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 457 455 ! 458 ! If time splitting, save data at first barotropic iteration 459 IF ( PRESENT(kit) ) THEN 460 IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data: 461 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1)) 462 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2)) 463 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3)) 464 465 ELSE ! Initialize arrays from slow varying open boundary data: 466 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 467 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 468 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 469 ENDIF 456 ! If time splitting, initialize arrays from slow varying open boundary data: 457 IF ( PRESENT(kit) ) THEN 458 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 459 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 460 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 470 461 ENDIF 471 462 ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
r5901 r6079 10 10 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 11 11 !!---------------------------------------------------------------------- 12 #if defined key_bdy && defined key_dynspg_flt12 #if defined key_bdy 13 13 !!---------------------------------------------------------------------- 14 !! 'key_bdy' AND unstructured open boundary conditions 15 !! 'key_dynspg_flt' filtered free surface 14 !! 'key_bdy' unstructured open boundary conditions 16 15 !!---------------------------------------------------------------------- 17 16 USE oce ! ocean dynamics and tracers … … 30 29 PRIVATE 31 30 32 PUBLIC bdy_vol ! routine called by dynspg_flt.h9031 PUBLIC bdy_vol 33 32 34 33 !! * Substitutions … … 45 44 !! *** ROUTINE bdyvol *** 46 45 !! 47 !! ** Purpose : This routine is called in dynspg_flt to control48 !! the volume of the system.A correction velocity is calculated46 !! ** Purpose : This routine controls the volume of the system. 47 !! A correction velocity is calculated 49 48 !! to correct the total transport through the unstructured OBC. 50 49 !! The total depth used is constant (H0) to be consistent with the -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90
r5901 r6079 18 18 #endif 19 19 USE dyncor_c1d ! Coriolis term (c1d case) (dyn_cor_1d ) 20 USE dynnxt _c1d! time-stepping (dyn_nxt routine)20 USE dynnxt ! time-stepping (dyn_nxt routine) 21 21 USE dyndmp ! U & V momentum damping (dyn_dmp routine) 22 22 USE restart ! restart … … 139 139 CALL dyn_cor_c1d( kstp ) ! vorticity term including Coriolis 140 140 CALL dyn_zdf ( kstp ) ! vertical diffusion 141 CALL dyn_nxt _c1d( kstp ) ! lateral velocity at next time step141 CALL dyn_nxt ( kstp ) ! lateral velocity at next time step 142 142 143 143 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r5620 r6079 14 14 USE dom_oce ! ocean space and time domain 15 15 USE phycst 16 USE dynspg_oce17 USE dynspg_ts18 16 USE daymod 19 17 USE tide_mod -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r5901 r6079 30 30 USE zdf_oce ! ocean vertical physics 31 31 USE ldftra ! lateral physics: eddy diffusivity coef. 32 USE sol_oce ! solver variables33 32 USE sbc_oce ! Surface boundary condition: ocean fields 34 33 USE sbc_ice ! Surface boundary condition: ice fields … … 46 45 USE iom 47 46 USE ioipsl 48 USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities49 47 50 48 #if defined key_lim2 … … 206 204 CALL iom_put( "sbu", z2d ) ! bottom i-current 207 205 ENDIF 208 #if defined key_dynspg_ts 209 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 210 #else 211 CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current 212 #endif 206 207 IF ( ln_dynspg_ts ) THEN 208 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 209 ELSE 210 CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current 211 ENDIF 213 212 214 213 CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current … … 223 222 CALL iom_put( "sbv", z2d ) ! bottom j-current 224 223 ENDIF 225 #if defined key_dynspg_ts 226 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current 227 #else 228 CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current 229 #endif 224 225 IF ( ln_dynspg_ts ) THEN 226 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current 227 ELSE 228 CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current 229 ENDIF 230 230 231 231 CALL iom_put( "woce", wn ) ! vertical velocity -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5901 r6079 45 45 LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers 46 46 47 !! Free surface parameters 48 !! ======================= 49 LOGICAL, PUBLIC :: ln_dynspg_exp !: Explicit free surface flag 50 LOGICAL, PUBLIC :: ln_dynspg_ts !: Split-Explicit free surface flag 51 47 52 !! Time splitting parameters 48 53 !! ========================= 49 54 LOGICAL, PUBLIC :: ln_bt_fw !: Forward integration of barotropic sub-stepping 50 55 LOGICAL, PUBLIC :: ln_bt_av !: Time averaging of barotropic variables 51 LOGICAL, PUBLIC :: ln_bt_ nn_auto!: Set number of barotropic iterations automatically56 LOGICAL, PUBLIC :: ln_bt_auto !: Set number of barotropic iterations automatically 52 57 INTEGER, PUBLIC :: nn_bt_flt !: Filter choice 53 58 INTEGER, PUBLIC :: nn_baro !: Number of barotropic iterations during one baroclinic step (rdt) 54 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_ nn_auto=T)59 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_auto=T) 55 60 56 61 !! Horizontal grid parameters for domhgr -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5901 r6079 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE lib_mpp 30 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient31 30 USE wrk_nemo ! Memory allocation 32 31 USE timing ! Timing -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5901 r6079 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 36 USE domvvl ! varying vertical mesh 37 USE dynspg_oce ! pressure gradient schemes38 USE dynspg_flt ! filtered free surface39 USE sol_oce ! ocean solver variables40 37 ! 41 38 USE in_out_manager ! I/O manager … … 133 130 ENDIF 134 131 ! 135 IF( lk_agrif ) THEN ! read free surface arrays in restart file136 IF( ln_rstart ) THEN137 IF( lk_dynspg_flt ) THEN ! read or initialize the following fields138 ! ! gcx, gcxb for agrif_opa_init139 IF( sol_oce_alloc() > 0 ) CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')140 CALL flt_rst( nit000, 'READ' )141 ENDIF142 ENDIF ! explicit case not coded yet with AGRIF143 ENDIF144 !145 132 ! 146 133 ! Initialize "now" and "before" barotropic velocities: … … 196 183 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)' 197 184 ! 198 #if ! defined key_istate_fixed199 185 DO jk = 1, jpk 200 186 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) & … … 202 188 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 203 189 END DO 204 #else205 tsn(:,:,:,jp_tem) = ztem * tmask(:,:,:)206 tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)207 #endif208 190 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 209 191 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 210 211 192 ! 212 193 END SUBROUTINE istate_t_s … … 451 432 !! p=integral [ rau*g dz ] 452 433 !!---------------------------------------------------------------------- 453 USE dynspg ! surface pressure gradient (dyn_spg routine)454 434 USE divhor ! hor. divergence (div_hor routine) 455 435 USE lbclnk ! ocean lateral boundary condition (or mpp link) 456 436 ! 457 437 INTEGER :: ji, jj, jk ! dummy loop indices 458 INTEGER :: indic ! ???459 438 REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars 460 439 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn … … 523 502 vb(:,:,:) = vn(:,:,:) 524 503 525 ! WARNING !!!!!526 ! after initializing u and v, we need to calculate the initial streamfunction bsf.527 ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).528 ! to do that, we call dyn_spg with a special trick:529 ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the530 ! right value assuming the velocities have been set up in one time step.531 ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)532 ! sets up s false trend to calculate the barotropic streamfunction.533 534 ua(:,:,:) = ub(:,:,:) / rdt535 va(:,:,:) = vb(:,:,:) / rdt536 537 ! calls dyn_spg. we assume euler time step, starting from rest.538 indic = 0539 CALL dyn_spg( nit000, indic ) ! surface pressure gradient540 !541 ! the new velocity is ua*rdt542 !543 CALL lbc_lnk( ua, 'U', -1. )544 CALL lbc_lnk( va, 'V', -1. )545 546 ub(:,:,:) = ua(:,:,:) * rdt547 vb(:,:,:) = va(:,:,:) * rdt548 ua(:,:,:) = 0.e0549 va(:,:,:) = 0.e0550 un(:,:,:) = ub(:,:,:)551 vn(:,:,:) = vb(:,:,:)552 504 ! 553 505 !!gm Check here call to div_hor should not be necessary -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5901 r6079 36 36 USE trd_oce ! trends: ocean variables 37 37 USE trddyn ! trend manager: dynamics 38 !jc USE zpshde ! partial step: hor. derivative (zps_hde routine) 38 39 ! 39 40 USE in_out_manager ! I/O manager … … 58 59 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 60 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag61 61 62 62 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) … … 132 132 !! 133 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf , ln_dynhpg_imp134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 135 135 !!---------------------------------------------------------------------- 136 136 ! … … 155 155 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 156 156 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 157 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp158 157 ENDIF 159 158 ! … … 293 292 ENDIF 294 293 294 ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 295 !jc CALL zps_hde ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 295 296 296 297 ! Local constant initialization … … 491 492 ! iniitialised to 0. zhpi zhpi 492 493 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 495 ! Partial steps: top & bottom before horizontal gradient 496 !jc CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi, & 497 !jc & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 498 !jc & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 493 499 494 500 !================================================================================== -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5901 r6079 19 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 20 20 !! 3.7 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 21 !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification 21 22 !!------------------------------------------------------------------------- 22 23 … … 28 29 USE sbc_oce ! Surface boundary condition: ocean fields 29 30 USE phycst ! physical constants 30 USE dynspg_oce ! type of surface pressure gradient31 31 USE dynadv ! dynamics: vector invariant versus flux form 32 32 USE domvvl ! variable volume … … 68 68 !! *** ROUTINE dyn_nxt *** 69 69 !! 70 !! ** Purpose : Compute the after horizontal velocity. Apply the boundary71 !! condition on the after velocity, achieve dthe time stepping70 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 71 !! condition on the after velocity, achieve the time stepping 72 72 !! by applying the Asselin filter on now fields and swapping 73 73 !! the fields. 74 74 !! 75 !! ** Method : * After velocity is compute using a leap-frog scheme: 76 !! (ua,va) = (ub,vb) + 2 rdt (ua,va) 77 !! Note that with flux form advection and variable volume layer 78 !! (lk_vvl=T), the leap-frog is applied on thickness weighted 79 !! velocity. 80 !! Note also that in filtered free surface (lk_dynspg_flt=T), 81 !! the time stepping has already been done in dynspg module 75 !! ** Method : * Ensure after velocities transport matches time splitting 76 !! estimate (ln_dynspg_ts=T) 82 77 !! 83 78 !! * Apply lateral boundary conditions on after velocity … … 101 96 INTEGER :: ji, jj, jk ! dummy loop indices 102 97 INTEGER :: iku, ikv ! local integers 103 #if ! defined key_dynspg_flt 104 REAL(wp) :: z2dt ! temporary scalar 105 #endif 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 107 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 108 100 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 112 104 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 113 105 ! 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 115 IF( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 106 IF( ln_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 107 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 108 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, zua, zva ) 116 109 ! 117 110 IF( kt == nit000 ) THEN … … 121 114 ENDIF 122 115 123 #if defined key_dynspg_flt 124 ! 125 ! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine 126 ! ------------- 127 128 ! Update after velocity on domain lateral boundaries (only local domain required) 129 ! -------------------------------------------------- 130 CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries 131 CALL lbc_lnk( va, 'V', -1. ) 132 ! 133 #else 134 135 # if defined key_dynspg_exp 136 ! Next velocity : Leap-frog time stepping 137 ! ------------- 138 z2dt = 2. * rdt ! Euler or leap-frog time step 139 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 140 ! 141 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 116 IF ( ln_dynspg_ts ) THEN 117 ! Ensure below that barotropic velocities match time splitting estimate 118 ! Compute actual transport and replace it with ts estimate at "after" time step 119 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 120 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 121 DO jk = 2, jpkm1 122 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 123 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 124 END DO 142 125 DO jk = 1, jpkm1 143 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 144 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 145 END DO 146 ELSE ! applied on thickness weighted velocity 147 DO jk = 1, jpkm1 148 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 149 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 150 & / fse3u_a(:,:,jk) * umask(:,:,jk) 151 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 152 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 153 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 154 END DO 155 ENDIF 156 # endif 157 158 # if defined key_dynspg_ts 159 !!gm IF ( lk_dynspg_ts ) THEN .... 160 ! Ensure below that barotropic velocities match time splitting estimate 161 ! Compute actual transport and replace it with ts estimate at "after" time step 162 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 163 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 164 DO jk = 2, jpkm1 165 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 167 END DO 168 DO jk = 1, jpkm1 169 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 171 END DO 172 173 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 174 ! Remove advective velocity from "now velocities" 175 ! prior to asselin filtering 176 ! In the forward case, this is done below after asselin filtering 177 ! so that asselin contribution is removed at the same time 178 DO jk = 1, jpkm1 179 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 180 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 181 END DO 182 ENDIF 183 !!gm ENDIF 184 # endif 126 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 127 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 128 END DO 129 130 IF ( .NOT.ln_bt_fw ) THEN 131 ! Remove advective velocity from "now velocities" 132 ! prior to asselin filtering 133 ! In the forward case, this is done below after asselin filtering 134 ! so that asselin contribution is removed at the same time 135 DO jk = 1, jpkm1 136 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 137 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 138 END DO 139 ENDIF 140 ENDIF 185 141 186 142 ! Update after velocity on domain lateral boundaries 187 143 ! -------------------------------------------------- 188 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries189 CALL lbc_lnk( va, 'V', -1. )190 !191 # if defined key_bdy192 ! !* BDY open boundaries193 IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )194 IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )195 196 !!$ Do we need a call to bdy_vol here??197 !198 # endif199 !200 144 # if defined key_agrif 201 145 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 202 146 # endif 203 #endif 204 147 ! 148 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries 149 CALL lbc_lnk( va, 'V', -1. ) 150 ! 151 # if defined key_bdy 152 ! !* BDY open boundaries 153 IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 154 IF( lk_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 155 156 !!$ Do we need a call to bdy_vol here?? 157 ! 158 # endif 159 ! 205 160 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 206 161 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step … … 259 214 ! (used as a now filtered scale factor until the swap) 260 215 ! ---------------------------------------------------- 261 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN216 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 262 217 ! No asselin filtering on thicknesses if forward time splitting 263 218 fse3t_b(:,:,:) = fse3t_n(:,:,:) 264 219 ELSE 265 220 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) … … 336 291 ENDIF 337 292 ! 338 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 339 294 ! Revert "before" velocities to time split estimate 340 295 ! Doing it here also means that asselin filter contribution is removed … … 400 355 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 401 356 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 402 ! 403 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 404 IF( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 357 ! 358 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 359 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 360 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 405 361 ! 406 362 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5901 r6079 6 6 !! History : 1.0 ! 2005-12 (C. Talandier, G. Madec, V. Garnier) Original code 7 7 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! dyn_spg : update the dynamics trend with the lateral diffusion 12 !! dyn_spg_ctl : initialization, namelist read, and parameters control 8 !! 3.7 ! 2015-11 (J. Chanut) Suppression of filtered free surface 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! dyn_spg : update the dynamics trend with surface pressure gradient 13 !! dyn_spg_init: initialization, namelist read, and parameters control 13 14 !!---------------------------------------------------------------------- 14 15 USE oce ! ocean dynamics and tracers variables … … 18 19 USE sbc_oce ! surface boundary condition: ocean 19 20 USE sbcapr ! surface boundary condition: atmospheric pressure 20 USE dynspg_oce ! surface pressure gradient variables21 21 USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine) 22 22 USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine) 23 USE dynspg_flt ! surface pressure gradient (dyn_spg_flt routine)24 USE dynadv ! dynamics: vector invariant versus flux form25 USE dynhpg, ONLY: ln_dynhpg_imp26 23 USE sbctide 27 24 USE updtide … … 32 29 USE in_out_manager ! I/O manager 33 30 USE lib_mpp ! MPP library 34 USE solver ! solver initialization35 31 USE wrk_nemo ! Memory Allocation 36 32 USE timing ! Timing … … 43 39 PUBLIC dyn_spg_init ! routine called by opa module 44 40 45 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from l k_dynspg_...41 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from ln_dynspg_... 46 42 47 43 !! * Substitutions … … 55 51 CONTAINS 56 52 57 SUBROUTINE dyn_spg( kt , kindic)53 SUBROUTINE dyn_spg( kt ) 58 54 !!---------------------------------------------------------------------- 59 55 !! *** ROUTINE dyn_spg *** … … 66 62 !!gm the after velocity, in the 2 other (ua,va) are still the trends 67 63 !! 68 !! ** Method : T hreeschemes:64 !! ** Method : Two schemes: 69 65 !! - explicit computation : the spg is evaluated at now 70 !! - filtered computation : the Roulet & madec (2000) technique is used71 66 !! - split-explicit computation: a time splitting technique is used 72 67 !! … … 80 75 ! 81 76 INTEGER, INTENT(in ) :: kt ! ocean time-step index 82 INTEGER, INTENT( out) :: kindic ! solver flag83 77 ! 84 78 INTEGER :: ji, jj, jk ! dummy loop indices … … 103 97 104 98 IF( ln_apr_dyn & ! atmos. pressure 105 .OR. ( .NOT.l k_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting)99 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting) 106 100 .OR. nn_ice_embd == 2 ) THEN ! embedded sea-ice 107 101 ! … … 113 107 END DO 114 108 ! 115 IF( ln_apr_dyn .AND. (.NOT. l k_dynspg_ts) ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==!109 IF( ln_apr_dyn .AND. (.NOT. ln_dynspg_ts) ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==! 116 110 zg_2 = grav * 0.5 117 111 DO jj = 2, jpjm1 ! gradient of Patm using inverse barometer ssh … … 126 120 ! 127 121 ! !== tide potential forcing term ==! 128 IF( .NOT.l k_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case122 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case 129 123 ! 130 124 CALL upd_tide( kt ) ! update tide potential … … 171 165 CASE ( 0 ) ; CALL dyn_spg_exp( kt ) ! explicit 172 166 CASE ( 1 ) ; CALL dyn_spg_ts ( kt ) ! time-splitting 173 CASE ( 2 ) ; CALL dyn_spg_flt( kt, kindic ) ! filtered174 167 ! 175 168 END SELECT 176 169 ! 177 170 IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics 178 SELECT CASE ( nspg ) 179 CASE ( 0, 1 ) 180 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 181 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 182 CASE( 2 ) 183 z2dt = 2. * rdt 184 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 185 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 186 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 187 END SELECT 171 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 188 173 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 189 174 ! … … 203 188 !! *** ROUTINE dyn_spg_init *** 204 189 !! 205 !! ** Purpose : Control the consistency between cppoptions for190 !! ** Purpose : Control the consistency between namelist options for 206 191 !! surface pressure gradient schemes 207 192 !!---------------------------------------------------------------------- 208 INTEGER :: ioptio 193 INTEGER :: ioptio, ios 194 ! 195 NAMELIST/namdyn_spg/ ln_dynspg_exp, ln_dynspg_ts, & 196 & ln_bt_fw, ln_bt_av, ln_bt_auto, & 197 & nn_baro, rn_bt_cmax, nn_bt_flt 209 198 !!---------------------------------------------------------------------- 210 199 ! 211 200 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_init') 212 201 ! 213 IF(lwp) THEN ! Control print 202 REWIND( numnam_ref ) ! Namelist namdyn_spg in reference namelist : Free surface 203 READ ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901) 204 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist', lwp ) 205 206 REWIND( numnam_cfg ) ! Namelist namdyn_spg in configuration namelist : Free surface 207 READ ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 ) 208 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist', lwp ) 209 IF(lwm) WRITE ( numond, namdyn_spg ) 210 ! 211 IF(lwp) THEN ! Namelist print 214 212 WRITE(numout,*) 215 213 WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme' 216 214 WRITE(numout,*) '~~~~~~~~~~~' 217 WRITE(numout,*) ' Explicit free surface lk_dynspg_exp = ', lk_dynspg_exp 218 WRITE(numout,*) ' Free surface with time splitting lk_dynspg_ts = ', lk_dynspg_ts 219 WRITE(numout,*) ' Filtered free surface cst volume lk_dynspg_flt = ', lk_dynspg_flt 220 ENDIF 221 222 IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 223 ! (do it now, to set nn_baro, used to allocate some arrays later on) 224 ! ! allocate dyn_spg arrays 225 IF( lk_dynspg_ts ) THEN 226 IF( dynspg_oce_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_oce arrays') 215 WRITE(numout,*) ' Explicit free surface ln_dynspg_exp = ', ln_dynspg_exp 216 WRITE(numout,*) ' Free surface with time splitting ln_dynspg_ts = ', ln_dynspg_ts 217 ENDIF 218 219 IF( ln_dynspg_ts ) THEN 220 CALL dyn_spg_ts_init( nit000 ) ! do it first, to set nn_baro, used to allocate some arrays later on 221 ! ! allocate dyn_spg arrays 227 222 IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays') 228 223 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) … … 231 226 ! ! Control of surface pressure gradient scheme options 232 227 ioptio = 0 233 IF(lk_dynspg_exp) ioptio = ioptio + 1 234 IF(lk_dynspg_ts ) ioptio = ioptio + 1 235 IF(lk_dynspg_flt) ioptio = ioptio + 1 228 IF(ln_dynspg_exp) ioptio = ioptio + 1 229 IF(ln_dynspg_ts ) ioptio = ioptio + 1 236 230 ! 237 231 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 238 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 239 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav ) & 240 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 241 ! 242 IF( lk_dynspg_exp) nspg = 0 243 IF( lk_dynspg_ts ) nspg = 1 244 IF( lk_dynspg_flt) nspg = 2 232 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 233 IF( ln_dynspg_ts .AND. ln_isfcav ) & 234 & CALL ctl_stop( ' dynspg_ts not tested with ice shelf cavity ' ) 235 ! 236 IF( ln_dynspg_exp) nspg = 0 237 IF( ln_dynspg_ts ) nspg = 1 245 238 ! 246 239 IF(lwp) THEN … … 248 241 IF( nspg == 0 ) WRITE(numout,*) ' explicit free surface' 249 242 IF( nspg == 1 ) WRITE(numout,*) ' free surface with time splitting scheme' 250 IF( nspg == 2 ) WRITE(numout,*) ' filtered free surface' 251 ENDIF 252 253 #if defined key_dynspg_flt 254 CALL solver_init( nit000 ) ! Elliptic solver initialisation 255 #endif 256 ! ! Control of hydrostatic pressure choice 257 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 243 ENDIF 258 244 ! 259 245 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_init') -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r5901 r6079 7 7 !! 3.2 ! 2009-06 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 8 8 !!---------------------------------------------------------------------- 9 #if defined key_dynspg_exp10 9 !!---------------------------------------------------------------------- 11 !! 'key_dynspg_exp'explicit free surface10 !! explicit free surface 12 11 !!---------------------------------------------------------------------- 13 12 !! dyn_spg_exp : update the momentum trend with the surface … … 31 30 PRIVATE 32 31 33 PUBLIC dyn_spg_exp ! routine called by step.F9032 PUBLIC dyn_spg_exp ! routine called by dyn_spg 34 33 35 34 !! * Substitutions … … 101 100 END SUBROUTINE dyn_spg_exp 102 101 103 #else104 !!----------------------------------------------------------------------105 !! Default case : Empty module No standart explicit free surface106 !!----------------------------------------------------------------------107 CONTAINS108 SUBROUTINE dyn_spg_exp( kt ) ! Empty routine109 WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt110 END SUBROUTINE dyn_spg_exp111 #endif112 113 102 !!====================================================================== 114 103 END MODULE dynspg_exp -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5901 r6079 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 13 14 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts15 15 !!---------------------------------------------------------------------- 16 !! 'key_dynspg_ts'split explicit free surface16 !! split explicit free surface 17 17 !!---------------------------------------------------------------------- 18 18 !! dyn_spg_ts : compute surface pressure gradient trend using a time- … … 23 23 USE sbc_oce ! surface boundary condition: ocean 24 24 USE sbcisf ! ice shelf variable (fwfisf) 25 USE dynspg_oce ! surface pressure gradient variables26 25 USE phycst ! physical constants 27 26 USE dynvor ! vorticity term 28 27 USE bdy_par ! for lk_bdy 29 USE bdytides ! open boundary condition data 28 USE bdytides ! open boundary condition data 30 29 USE bdydyn2d ! open boundary conditions on barotropic variables 31 30 USE sbctide ! tides … … 68 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 69 68 70 ! Arrays below are saved to allow testing of the "no time averaging" option71 ! If this option is not retained, these could be replaced by temporary arrays72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays73 ubb_e, ub_e, &74 vbb_e, vb_e75 76 69 !! * Substitutions 77 70 # include "domzgr_substitute.h90" … … 88 81 !! *** routine dyn_spg_ts_alloc *** 89 82 !!---------------------------------------------------------------------- 90 INTEGER :: ierr( 3)83 INTEGER :: ierr(4) 91 84 !!---------------------------------------------------------------------- 92 85 ierr(:) = 0 93 86 94 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 95 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 96 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 87 ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 88 & ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), & 89 & va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), & 90 & hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT= ierr(1) ) 97 91 98 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) … … 101 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 102 96 97 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 98 #if defined key_agrif 99 & ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , & 100 #endif 101 & STAT= ierr(4)) 102 103 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 104 105 105 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn spg_oce_alloc: failed to allocate arrays')106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 107 107 ! 108 108 END FUNCTION dyn_spg_ts_alloc … … 146 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 147 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - -149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 150 REAL(wp) :: zu_spg, zv_spg ! - -151 REAL(wp) :: zhura, zhvra 152 REAL(wp) :: za0, za1, za2, za3 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: z un_e, zvn_e, zsshp2_e148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 150 REAL(wp) :: zu_spg, zv_spg ! - - 151 REAL(wp) :: zhura, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 ! - - 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 155 155 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 156 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv156 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 157 157 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 158 158 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a … … 164 164 ! !* Allocate temporary arrays 165 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)167 CALL wrk_alloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc)166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) … … 188 188 ! 189 189 ! time offset in steps for bdy data update 190 IF (.NOT.ln_bt_fw) THEN ; noffset=- 2*nn_baro ; ELSE ; noffset = 0 ; ENDIF190 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 191 191 ! 192 192 IF( kt == nit000 ) THEN !* initialisation … … 485 485 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 486 486 sshn_e(:,:) = sshn (:,:) 487 zun_e(:,:) = un_b (:,:)488 zvn_e(:,:) = vn_b (:,:)487 un_e (:,:) = un_b (:,:) 488 vn_e (:,:) = vn_b (:,:) 489 489 ! 490 490 hu_e (:,:) = hu (:,:) … … 494 494 ELSE ! CENTRED integration: start from BEFORE fields 495 495 sshn_e(:,:) = sshb (:,:) 496 zun_e(:,:) = ub_b (:,:)497 zvn_e(:,:) = vb_b (:,:)496 un_e (:,:) = ub_b (:,:) 497 vn_e (:,:) = vb_b (:,:) 498 498 ! 499 499 hu_e (:,:) = hu_b (:,:) … … 509 509 va_b (:,:) = 0._wp 510 510 ssha (:,:) = 0._wp ! Sum for after averaged sea level 511 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop512 zv_sum(:,:) = 0._wp511 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 512 vn_adv(:,:) = 0._wp 513 513 ! ! ==================== ! 514 514 DO jn = 1, icycle ! sub-time-step loop ! … … 518 518 ! Update only tidal forcing at open boundaries 519 519 #if defined key_tide 520 IF ( lk_bdy .AND. lk_tide ) 521 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )520 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 521 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 522 522 #endif 523 523 ! … … 534 534 535 535 ! Extrapolate barotropic velocities at step jit+0.5: 536 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)537 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)536 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 537 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 538 538 539 539 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) … … 597 597 ! Sum over sub-time-steps to compute advective velocities 598 598 za2 = wgtbtp2(jn) 599 zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)600 zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)599 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 600 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 601 601 ! 602 602 ! Set next sea level: … … 733 733 ! 734 734 ! Add bottom stresses: 735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 737 ! 738 738 ! Surface pressure trend: … … 751 751 DO jj = 2, jpjm1 752 752 DO ji = fs_2, fs_jpim1 ! vector opt. 753 ua_e(ji,jj) = ( zun_e(ji,jj) &753 ua_e(ji,jj) = ( un_e(ji,jj) & 754 754 & + rdtbt * ( zwx(ji,jj) & 755 755 & + zu_trd(ji,jj) & … … 757 757 & ) * umask(ji,jj,1) 758 758 759 va_e(ji,jj) = ( zvn_e(ji,jj) &759 va_e(ji,jj) = ( vn_e(ji,jj) & 760 760 & + rdtbt * ( zwy(ji,jj) & 761 761 & + zv_trd(ji,jj) & … … 772 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 773 773 774 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &774 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 775 775 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 776 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & … … 778 778 & ) * zhura 779 779 780 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &780 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 781 781 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 782 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & … … 802 802 #if defined key_bdy 803 803 ! open boundaries 804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 805 805 #endif 806 806 #if defined key_agrif … … 810 810 ! ! ---- 811 811 ubb_e (:,:) = ub_e (:,:) 812 ub_e (:,:) = zun_e(:,:)813 zun_e(:,:) = ua_e (:,:)812 ub_e (:,:) = un_e (:,:) 813 un_e (:,:) = ua_e (:,:) 814 814 ! 815 815 vbb_e (:,:) = vb_e (:,:) 816 vb_e (:,:) = zvn_e(:,:)817 zvn_e(:,:) = va_e (:,:)816 vb_e (:,:) = vn_e (:,:) 817 vn_e (:,:) = va_e (:,:) 818 818 ! 819 819 sshbb_e(:,:) = sshb_e(:,:) … … 840 840 ! ----------------------------------------------------------------------------- 841 841 ! 842 ! At this stage ssha holds a time averaged value 843 ! ! Sea Surface Height at u-,v- and f-points 844 IF( lk_vvl ) THEN ! (required only in key_vvl case) 842 ! Set advection velocity correction: 843 zwx(:,:) = un_adv(:,:) 844 zwy(:,:) = vn_adv(:,:) 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 846 un_adv(:,:) = zwx(:,:)*hur(:,:) 847 vn_adv(:,:) = zwy(:,:)*hvr(:,:) 848 ELSE 849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 851 END IF 852 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 854 ub2_b(:,:) = zwx(:,:) 855 vb2_b(:,:) = zwy(:,:) 856 ENDIF 857 ! 858 ! Update barotropic trend: 859 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 860 DO jk=1,jpkm1 861 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 862 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 863 END DO 864 ELSE 865 ! At this stage, ssha has been corrected: compute new depths at velocity points 845 866 DO jj = 1, jpjm1 846 867 DO ji = 1, jpim1 ! NO Vector Opt. … … 854 875 END DO 855 876 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 856 ENDIF 857 ! 858 ! Set advection velocity correction: 859 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 860 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 861 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 862 ELSE 863 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 864 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 865 END IF 866 867 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 868 ub2_b(:,:) = zu_sum(:,:) 869 vb2_b(:,:) = zv_sum(:,:) 870 ENDIF 871 ! 872 ! Update barotropic trend: 873 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 874 DO jk=1,jpkm1 875 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 876 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 877 END DO 878 ELSE 877 ! 879 878 DO jk=1,jpkm1 880 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 915 914 ! 916 915 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 917 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)918 CALL wrk_dealloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc )916 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 917 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 919 918 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 920 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) … … 1064 1063 ! 1065 1064 INTEGER :: ji ,jj 1066 INTEGER :: ios ! Local integer output status for namelist read1067 1065 REAL(wp) :: zxr2, zyr2, zcmax 1068 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1069 1067 !! 1070 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &1071 & nn_baro, rn_bt_cmax, nn_bt_flt1072 1068 !!---------------------------------------------------------------------- 1073 1069 ! 1074 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1075 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1076 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1077 1078 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1079 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1080 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1081 IF(lwm) WRITE ( numond, namsplit ) 1082 ! 1083 ! ! Max courant number for ext. grav. waves 1070 ! Max courant number for ext. grav. waves 1084 1071 ! 1085 1072 CALL wrk_alloc( jpi, jpj, zcu ) 1086 1073 ! 1087 IF (lk_vvl) THEN 1088 DO jj = 1, jpj 1089 DO ji =1, jpi 1090 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1091 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1092 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1093 END DO 1094 END DO 1095 ELSE 1096 DO jj = 1, jpj 1097 DO ji =1, jpi 1098 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1099 zyr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1100 zcu(ji,jj) = SQRT( grav * ht(ji,jj) * (zxr2 + zyr2) ) 1101 END DO 1102 END DO 1103 ENDIF 1104 1074 DO jj = 1, jpj 1075 DO ji =1, jpi 1076 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1077 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1078 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1079 END DO 1080 END DO 1081 ! 1105 1082 zcmax = MAXVAL( zcu(:,:) ) 1106 1083 IF( lk_mpp ) CALL mpp_max( zcmax ) 1107 1084 1108 1085 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1109 IF (ln_bt_ nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1086 IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1110 1087 1111 1088 rdtbt = rdt / REAL( nn_baro , wp ) … … 1115 1092 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1116 1093 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1117 IF( ln_bt_ nn_auto ) THEN1118 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1094 IF( ln_bt_auto ) THEN 1095 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1119 1096 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1120 1097 ELSE 1121 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1098 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1122 1099 ENDIF 1123 1100 … … 1164 1141 END SUBROUTINE dyn_spg_ts_init 1165 1142 1166 #else1167 !!---------------------------------------------------------------------------1168 !! Default case : Empty module No split explicit free surface1169 !!---------------------------------------------------------------------------1170 CONTAINS1171 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1172 dyn_spg_ts_alloc = 01173 END FUNCTION dyn_spg_ts_alloc1174 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1175 INTEGER, INTENT(in) :: kt1176 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1177 END SUBROUTINE dyn_spg_ts1178 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1179 INTEGER , INTENT(in) :: kt ! ocean time-step1180 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1181 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1182 END SUBROUTINE ts_rst1183 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1184 INTEGER , INTENT(in) :: kt ! ocean time-step1185 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1186 END SUBROUTINE dyn_spg_ts_init1187 #endif1188 1189 1143 !!====================================================================== 1190 1144 END MODULE dynspg_ts -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5901 r6079 32 32 USE trd_oce ! trends: ocean variables 33 33 USE trddyn ! trend manager: dynamics 34 USE c1d ! 1D vertical configuration 34 35 ! 35 36 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 547 548 END SELECT 548 549 ! 550 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 551 DO jj = 1, jpjm1 552 DO ji = 1, fs_jpim1 ! vector opt. 553 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 554 END DO 555 END DO 556 ENDIF 557 ! 549 558 CALL lbc_lnk( zwz, 'F', 1. ) 550 !551 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!552 DO jj = 1, jpjm1553 DO ji = 1, fs_jpim1 ! vector opt.554 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)555 END DO556 END DO557 ENDIF558 559 ! 559 560 ! !== horizontal fluxes ==! … … 662 663 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 663 664 ! 664 IF( ioptio /= 1) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )665 IF( ( ioptio /= 1).AND.( .NOT.lk_c1d ) ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 665 666 ! 666 667 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r3625 r6079 8 8 !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 !! 3.7 ! 2015-11 (J. Chanut) output velocities instead of trends 10 11 !!---------------------------------------------------------------------- 11 12 … … 18 19 USE phycst ! physical constants 19 20 USE zdf_oce ! ocean vertical physics 21 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 20 22 USE sbc_oce ! surface boundary condition: ocean 21 23 USE lib_mpp ! MPP library … … 118 120 ! 119 121 END DO ! End of time splitting 122 123 ! Time step momentum here to be compliant with what is done in the implicit case 124 ! 125 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 126 DO jk = 1, jpkm1 127 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 128 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 129 END DO 130 ELSE ! applied on thickness weighted velocity 131 DO jk = 1, jpkm1 132 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 133 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 134 & / fse3u_a(:,:,jk) * umask(:,:,jk) 135 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 136 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 137 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 138 END DO 139 ENDIF 120 140 ! 121 141 CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r5901 r6079 25 25 USE wrk_nemo ! Memory Allocation 26 26 USE timing ! Timing 27 USE dynadv ! dynamics: vector invariant versus flux form 28 USE dynspg_oce, ONLY: lk_dynspg_ts 27 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 29 28 30 29 IMPLICIT NONE … … 87 86 ENDIF 88 87 89 ! 0. Local constant initialization 90 ! -------------------------------- 91 z1_p2dt = 1._wp / p2dt ! inverse of the timestep 92 93 ! 1. Apply semi-implicit bottom friction 88 ! 1. Time step dynamics 89 ! --------------------- 90 91 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 92 DO jk = 1, jpkm1 93 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 94 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 95 END DO 96 ELSE ! applied on thickness weighted velocity 97 DO jk = 1, jpkm1 98 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 99 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 100 & / fse3u_a(:,:,jk) * umask(:,:,jk) 101 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 102 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 103 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 104 END DO 105 ENDIF 106 107 ! 2. Apply semi-implicit bottom friction 94 108 ! -------------------------------------- 95 109 ! Only needed for semi-implicit bottom friction setup. The explicit … … 97 111 ! column vector of the tri-diagonal matrix equation 98 112 ! 99 100 113 IF( ln_bfrimp ) THEN 101 114 DO jj = 2, jpjm1 … … 119 132 ENDIF 120 133 121 #if defined key_dynspg_ts 122 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 123 DO jk = 1, jpkm1 124 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 125 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 126 END DO 127 ELSE ! applied on thickness weighted velocity 128 DO jk = 1, jpkm1 129 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 130 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 131 & / fse3u_a(:,:,jk) * umask(:,:,jk) 132 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 133 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 134 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 135 END DO 136 ENDIF 137 138 IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 134 ! With split-explicit free surface, barotropic stress is treated explicitly 135 ! Update velocities at the bottom. 136 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 137 ! not lead to the effective stress seen over the whole barotropic loop. 138 IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN 139 139 ! remove barotropic velocities: 140 140 DO jk = 1, jpkm1 … … 166 166 END IF 167 167 ENDIF 168 #endif 169 170 ! 2. Vertical diffusion on u 168 169 ! 3. Vertical diffusion on u 171 170 ! --------------------------- 172 171 ! Matrix and second member construction … … 219 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 220 219 DO ji = fs_2, fs_jpim1 ! vector opt. 221 #if defined key_dynspg_ts222 220 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 223 221 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 224 222 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 225 #else226 ua(ji,jj,1) = ub(ji,jj,1) &227 & + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &228 & / ( fse3u(ji,jj,1) * rau0 ) * umask(ji,jj,1) )229 #endif230 223 END DO 231 224 END DO … … 233 226 DO jj = 2, jpjm1 234 227 DO ji = fs_2, fs_jpim1 235 #if defined key_dynspg_ts236 228 zrhs = ua(ji,jj,jk) ! zrhs=right hand side 237 #else238 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)239 #endif240 229 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 241 230 END DO … … 256 245 END DO 257 246 258 #if ! defined key_dynspg_ts 259 ! Normalization to obtain the general momentum trend ua 260 DO jk = 1, jpkm1 261 DO jj = 2, jpjm1 262 DO ji = fs_2, fs_jpim1 ! vector opt. 263 ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 264 END DO 265 END DO 266 END DO 267 #endif 268 269 ! 3. Vertical diffusion on v 247 ! 4. Vertical diffusion on v 270 248 ! --------------------------- 271 249 ! Matrix and second member construction … … 317 295 ! 318 296 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 319 DO ji = fs_2, fs_jpim1 ! vector opt. 320 #if defined key_dynspg_ts 297 DO ji = fs_2, fs_jpim1 ! vector opt. 321 298 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1) 322 299 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 323 300 & / ( ze3va * rau0 ) 324 #else325 va(ji,jj,1) = vb(ji,jj,1) &326 & + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &327 & / ( fse3v(ji,jj,1) * rau0 ) )328 #endif329 301 END DO 330 302 END DO … … 332 304 DO jj = 2, jpjm1 333 305 DO ji = fs_2, fs_jpim1 ! vector opt. 334 #if defined key_dynspg_ts335 306 zrhs = va(ji,jj,jk) ! zrhs=right hand side 336 #else337 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)338 #endif339 307 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 340 308 END DO … … 355 323 END DO 356 324 357 ! Normalization to obtain the general momentum trend va358 #if ! defined key_dynspg_ts359 DO jk = 1, jpkm1360 DO jj = 2, jpjm1361 DO ji = fs_2, fs_jpim1 ! vector opt.362 va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt363 END DO364 END DO365 END DO366 #endif367 368 325 ! J. Chanut: Lines below are useless ? 369 !! restore bottom layer avmu(v)326 !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces) 370 327 IF( ln_bfrimp ) THEN 371 328 DO jj = 2, jpjm1 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5901 r6079 106 106 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 107 107 108 #if ! defined key_dynspg_ts 109 ! These lines are not necessary with time splitting since110 ! boundary condition on sea level is set during ts loop108 IF ( .NOT.ln_dynspg_ts ) THEN 109 ! These lines are not necessary with time splitting since 110 ! boundary condition on sea level is set during ts loop 111 111 # if defined key_agrif 112 CALL agrif_ssh( kt )112 CALL agrif_ssh( kt ) 113 113 # endif 114 114 # if defined key_bdy 115 IF( lk_bdy ) THEN116 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary117 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries118 ENDIF115 IF( lk_bdy ) THEN 116 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 117 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 118 ENDIF 119 119 # endif 120 #endif 120 ENDIF 121 121 122 122 #if defined key_asminc … … 250 250 ENDIF 251 251 252 # if defined key_dynspg_ts 253 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 254 # else 255 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 256 #endif 252 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 253 !** Euler time-stepping: no filter 257 254 sshb(:,:) = sshn(:,:) ! before <-- now 258 255 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r5891 r6079 32 32 PUBLIC fld_map ! routine called by tides_init 33 33 PUBLIC fld_read, fld_fill ! called by sbc... modules 34 PUBLIC fld_clopn 34 35 35 36 TYPE, PUBLIC :: FLD_N !: Namelist field informations … … 297 298 ztinta = REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp ) 298 299 ztintb = 1. - ztinta 299 !CDIR COLLAPSE300 300 sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2) 301 301 ELSE ! nothing to do... … … 1214 1214 imonth = kmonth 1215 1215 iday = kday 1216 IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week 1217 isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 1218 llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month 1219 llprevyr = llprevmth .AND. nmonth == 1 1220 iyear = nyear - COUNT((/llprevyr /)) 1221 imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)) 1222 iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday) 1223 ENDIF 1216 1224 ELSE ! use current day values 1217 1225 IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week … … 1313 1321 & CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 1314 1322 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 1315 sdf(jf)%igrd = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init1316 sdf(jf)%ibdy = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init1317 1323 END DO 1318 1324 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
r5901 r6079 12 12 USE dom_oce ! ocean space and time domain 13 13 USE sbc_oce ! surface boundary condition 14 USE dynspg_oce ! surface pressure gradient variables15 14 USE phycst ! physical constants 16 15 USE fldread ! read input fields … … 112 111 IF(lwp) WRITE(numout,*) ' Inverse barometer added to OBC ssh data' 113 112 ENDIF 114 IF( ( ln_apr_obc ) .AND. .NOT. lk_dynspg_ts ) & 115 CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY possible with time-splitting' ) 113 !jc: stop below should rather be a warning 116 114 IF( ( ln_apr_obc ) .AND. .NOT. ln_apr_dyn ) & 117 115 CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY requires ln_apr_dyn=T' ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
r5620 r6079 10 10 USE phycst 11 11 USE daymod 12 USE dynspg_oce13 12 USE tideini 14 13 ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
r5620 r6079 10 10 USE phycst 11 11 USE daymod 12 USE dynspg_oce13 12 USE tide_mod 14 13 ! … … 96 95 & CALL ctl_stop('rdttideramp must be positive') 97 96 ! 98 IF( .NOT. lk_dynspg_ts ) CALL ctl_warn( 'sbc_tide : use of time splitting is recommended' )99 97 ! 100 98 ALLOCATE( ntide(nb_harmo) ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90
r5620 r6079 31 31 CONTAINS 32 32 33 SUBROUTINE upd_tide( kt, kit, kbaro, koffset )33 SUBROUTINE upd_tide( kt, kit, time_offset ) 34 34 !!---------------------------------------------------------------------- 35 35 !! *** ROUTINE upd_tide *** … … 42 42 !!---------------------------------------------------------------------- 43 43 INTEGER, INTENT(in) :: kt ! ocean time-step index 44 INTEGER, INTENT(in), OPTIONAL :: kit ! external mode sub-time-step index (lk_dynspg_ts=T only)45 INTEGER, INTENT(in), OPTIONAL :: kbaro ! number of sub-time-step (lk_dynspg_ts=T only)46 INTEGER, INTENT(in), OPTIONAL :: koffset ! time offset in number47 ! of sub-time-steps (lk_dynspg_ts=T only)44 INTEGER, INTENT(in), OPTIONAL :: kit ! external mode sub-time-step index (lk_dynspg_ts=T) 45 INTEGER, INTENT(in), OPTIONAL :: time_offset ! time offset in number 46 ! of internal steps (lk_dynspg_ts=F) 47 ! of external steps (lk_dynspg_ts=T) 48 48 ! 49 49 INTEGER :: joffset ! local integer … … 57 57 ! 58 58 joffset = 0 59 IF( PRESENT( koffset ) ) joffset = koffset59 IF( PRESENT( time_offset ) ) joffset = time_offset 60 60 ! 61 IF( PRESENT( kit ) .AND. PRESENT( kbaro )) THEN62 zt = zt + ( kit + 0.5_wp * ( joffset - 1 ) ) * rdt / REAL( kbaro, wp )61 IF( PRESENT( kit ) ) THEN 62 zt = zt + ( kit + joffset - 1 ) * rdt / REAL( nn_baro, wp ) 63 63 ELSE 64 64 zt = zt + joffset * rdt … … 74 74 IF( ln_tide_ramp ) THEN ! linear increase if asked 75 75 zt = ( kt - nit000 ) * rdt 76 IF( PRESENT( kit ) .AND. PRESENT( kbaro ) ) zt = zt + kit * rdt / REAL( kbaro, wp )76 IF( PRESENT( kit ) ) zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp ) 77 77 zramp = MIN( MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp ) 78 78 pot_astro(:,:) = zramp * pot_astro(:,:) … … 86 86 !!---------------------------------------------------------------------- 87 87 CONTAINS 88 SUBROUTINE upd_tide( kt, kit, kbaro, koffset )! Empty routine88 SUBROUTINE upd_tide( kt, kit, time_offset ) ! Empty routine 89 89 INTEGER, INTENT(in) :: kt ! integer arg, dummy routine 90 90 INTEGER, INTENT(in), OPTIONAL :: kit ! optional arg, dummy routine 91 INTEGER, INTENT(in), OPTIONAL :: kbaro ! optional arg, dummy routine 92 INTEGER, INTENT(in), OPTIONAL :: koffset ! optional arg, dummy routine 91 INTEGER, INTENT(in), OPTIONAL :: time_offset ! optional arg, dummy routine 93 92 WRITE(*,*) 'upd_tide: You should not have seen this print! error?', kt 94 93 END SUBROUTINE upd_tide -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r5901 r6079 26 26 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. 27 27 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 28 USE c1d ! 1D vertical configuration 28 29 ! 29 30 USE in_out_manager ! I/O manager … … 214 215 CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' ) 215 216 ENDIF 216 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 217 IF( (ioptio /= 1).AND. (.NOT. lk_c1d ) ) & 218 CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 217 219 ! 218 220 IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 ) & ! Centered -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
r5901 r6079 19 19 USE trd_oce ! trends: ocean variables 20 20 USE trdtra ! tracers trends 21 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient22 21 USE diaptr ! poleward transport diagnostics 23 22 ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
r5901 r6079 21 21 USE trd_oce ! trends: ocean variables 22 22 USE trdtra ! tracers trends manager 23 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient24 23 USE sbcrnf ! river runoffs 25 24 USE diaptr ! poleward transport diagnostics -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r5901 r6079 20 20 USE trd_oce ! trends: ocean variables 21 21 USE trdtra ! trends manager: tracers 22 USE dynspg_oce ! surface pressure gradient variables23 22 USE diaptr ! poleward transport diagnostics 24 23 ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r5901 r6079 18 18 USE traadv_fct ! acces to routine interp_4th_cpt 19 19 USE trdtra ! trends manager: tracers 20 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient21 20 USE diaptr ! poleward transport diagnostics 22 21 ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r5901 r6079 31 31 USE zdf_oce ! ocean vertical mixing 32 32 USE domvvl ! variable volume 33 USE dynspg_oce ! surface pressure gradient variables34 USE dynhpg ! hydrostatic pressure gradient35 33 USE trd_oce ! trends: ocean variables 36 34 USE trdtra ! trends manager: tracers … … 58 56 PUBLIC tra_nxt_vvl ! to be used in trcnxt 59 57 60 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg61 58 62 59 !! * Substitutions … … 90 87 !! 91 88 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 92 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)89 !! 93 90 !!---------------------------------------------------------------------- 94 91 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 105 102 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 106 103 IF(lwp) WRITE(numout,*) '~~~~~~~' 107 !108 rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp) ! Brown & Campana parameter for semi-implicit hpg109 104 ENDIF 110 105 111 106 ! Update after tracer on domain lateral boundaries 112 107 ! 108 ! 113 109 #if defined key_agrif 114 110 CALL Agrif_tra ! AGRIF zoom boundaries … … 181 177 !! 182 178 !! ** Method : - Apply a Asselin time filter on now fields. 183 !! - save in (ta,sa) an average over the three time levels184 !! which will be used to compute rdn and thus the semi-implicit185 !! hydrostatic pressure gradient (ln_dynhpg_imp = T)186 179 !! - swap tracer fields to prepare the next time_step. 187 !! This can be summurized for tempearture as:188 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T189 !! ztm = 0 otherwise190 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp)191 !! tb = tn + atfp*[ tb - 2 tn + ta ]192 !! tn = ta193 !! ta = ztm (NB: reset to 0 after eos_bn2 call)194 180 !! 195 181 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 196 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)182 !! 197 183 !!---------------------------------------------------------------------- 198 184 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 205 191 ! 206 192 INTEGER :: ji, jj, jk, jn ! dummy loop indices 207 LOGICAL :: ll_tra_hpg ! local logical208 193 REAL(wp) :: ztn, ztd ! local scalars 209 194 !!---------------------------------------------------------------------- … … 215 200 ENDIF 216 201 ! 217 IF( cdtype == 'TRA' ) THEN ; ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg218 ELSE ; ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg219 ENDIF220 202 ! 221 203 DO jn = 1, kjpt … … 230 212 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 231 213 ! 232 IF( ll_tra_hpg ) pta(ji,jj,jk,jn) = ztn + rbcp * ztd ! pta <-- Brown & Campana average233 214 END DO 234 215 END DO … … 248 229 !! 249 230 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 250 !! - save in (ta,sa) a thickness weighted average over the three251 !! time levels which will be used to compute rdn and thus the semi-252 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)253 231 !! - swap tracer fields to prepare the next time_step. 254 232 !! This can be summurized for tempearture as: 255 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T256 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] )257 !! ztm = 0 otherwise258 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )259 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] )260 !! tn = ta261 !! ta = zt (NB: reset to 0 after eos_bn2 call)262 233 !! 263 234 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 264 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)235 !! 265 236 !!---------------------------------------------------------------------- 266 237 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 276 247 277 248 !! 278 LOGICAL :: ll_tra _hpg, ll_traqsr, ll_rnf, ll_isf ! local logical249 LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical 279 250 INTEGER :: ji, jj, jk, jn ! dummy loop indices 280 251 REAL(wp) :: zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar … … 289 260 ! 290 261 IF( cdtype == 'TRA' ) THEN 291 ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg292 262 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 293 263 ll_rnf = ln_rnf ! active tracers case and river runoffs … … 298 268 END IF 299 269 ELSE 300 ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg301 270 ll_traqsr = .FALSE. ! active tracers case and NO solar penetration 302 271 ll_rnf = .FALSE. ! passive tracers or NO river runoffs … … 356 325 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 357 326 ! 358 IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only)359 ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d )360 pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average361 ENDIF362 327 END DO 363 328 END DO -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r5901 r6079 18 18 USE zdf_oce ! ocean vertical physics variables 19 19 USE sbc_oce ! surface boundary condition: ocean 20 USE dynspg_oce21 20 ! 22 21 USE ldftra ! lateral diffusion: eddy diffusivity -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90
r5620 r6079 71 71 INTEGER, PUBLIC, PARAMETER :: jpdyn_bfri = 12 !: implicit bottom friction (ln_bfrimp=.TRUE.) 72 72 INTEGER, PUBLIC, PARAMETER :: jpdyn_ken = 13 !: use for calculation of KE 73 INTEGER, PUBLIC, PARAMETER :: jpdyn_spgflt = 14 !: filter contribution to surface pressure gradient (spg_flt)74 INTEGER, PUBLIC, PARAMETER :: jpdyn_spgexp = 15 !: explicit contribution to surface pressure gradient (spg_flt)75 73 ! 76 74 !!---------------------------------------------------------------------- -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90
r5620 r6079 113 113 CASE( jpdyn_spg ) ; CALL iom_put( "utrd_spg", putrd ) ! surface pressure gradient 114 114 CALL iom_put( "vtrd_spg", pvtrd ) 115 CASE( jpdyn_spgexp ); CALL iom_put( "utrd_spgexp", putrd ) ! surface pressure gradient (explicit)116 CALL iom_put( "vtrd_spgexp", pvtrd )117 CASE( jpdyn_spgflt ); CALL iom_put( "utrd_spgflt", putrd ) ! surface pressure gradient (filtered)118 CALL iom_put( "vtrd_spgflt", pvtrd )119 115 CASE( jpdyn_pvo ) ; CALL iom_put( "utrd_pvo", putrd ) ! planetary vorticity 120 116 CALL iom_put( "vtrd_pvo", pvtrd ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90
r5901 r6079 120 120 CASE( jpdyn_hpg ) ; CALL iom_put( "ketrd_hpg", zke ) ! hydrostatic pressure gradient 121 121 CASE( jpdyn_spg ) ; CALL iom_put( "ketrd_spg", zke ) ! surface pressure gradient 122 CASE( jpdyn_spgexp ); CALL iom_put( "ketrd_spgexp", zke ) ! surface pressure gradient (explicit)123 CASE( jpdyn_spgflt ); CALL iom_put( "ketrd_spgflt", zke ) ! surface pressure gradient (filter)124 122 CASE( jpdyn_pvo ) ; CALL iom_put( "ketrd_pvo", zke ) ! planetary vorticity 125 123 CASE( jpdyn_rvo ) ; CALL iom_put( "ketrd_rvo", zke ) ! relative vorticity (or metric term) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/oce.F90
r5901 r6079 21 21 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ub , un , ua !: i-horizontal velocity [m/s] 22 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vb , vn , va !: j-horizontal velocity [m/s] 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ua_sv, va_sv !: Saved trends (time spliting) [m/s2]24 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wn !: vertical velocity [m/s] 25 24 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn !: horizontal divergence [s-1] … … 38 37 ! 39 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: spgu, spgv !: horizontal surface pressure gradient 39 40 !! Specific to split explicit free surface (allocated in dynspg_ts module): 41 ! 42 !! Time filtered arrays at baroclinic time step: 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv, vn_adv ! Advection vel. at "now" barocl. step 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b , vb2_b ! Half step fluxes (ln_bt_fw=T) 45 #if defined key_agrif 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b ! Half step time integrated fluxes 47 #endif 48 49 !! Arrays at barotropic time step: ! bef before ! before ! now ! after ! 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubb_e , ub_e , un_e , ua_e 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vbb_e , vb_e , vn_e , va_e 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e , sshb_e , sshn_e , ssha_e 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_e 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_e 55 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur_e 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hvr_e 40 57 41 58 !! interpolated gradient (only used in zps case) … … 85 102 ! 86 103 ALLOCATE( ub (jpi,jpj,jpk) , un (jpi,jpj,jpk) , ua(jpi,jpj,jpk) , & 87 & vb (jpi,jpj,jpk) , vn (jpi,jpj,jpk) , va(jpi,jpj,jpk) , & 88 & ua_sv(jpi,jpj,jpk) , va_sv(jpi,jpj,jpk) , & 104 & vb (jpi,jpj,jpk) , vn (jpi,jpj,jpk) , va(jpi,jpj,jpk) , & 89 105 & wn (jpi,jpj,jpk) , hdivn(jpi,jpj,jpk) , & 90 106 & tsb (jpi,jpj,jpk,jpts) , tsn (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) , & -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/step.F90
r5901 r6079 28 28 !! 3.7 ! 2014-10 (G. Madec) LDF simplication 29 29 !! - ! 2014-12 (G. Madec) remove KPP scheme 30 !! - ! 2015-11 (J. Chanut) free surface simplification 30 31 !!---------------------------------------------------------------------- 31 32 … … 179 180 180 181 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 181 ! Ocean dynamics : hdiv, ssh, e3, wn 182 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 182 ! Ocean dynamics : hdiv, ssh, e3, u, v, w 183 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 184 183 185 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 184 186 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 185 187 CALL wzv ( kstp ) ! now cross-level velocity 186 188 187 IF( lk_dynspg_ts ) THEN188 ! In case the time splitting case, update almost all momentum trends here:189 ! Note that the computation of vertical velocity above, hence "after" sea level190 ! is necessary to compute momentum advection for the rhs of barotropic loop:191 189 !!gm : why also here ???? 192 IF(ln_sto_eos )CALL sto_pts( tsn ) ! Random T/S fluctuations190 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 193 191 !!gm 194 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 195 192 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 193 194 !!jc: fs simplification 195 !!jc: lines below are useless if lk_vvl=T. Keep them here (which maintains a bug if lk_vvl=F and ln_zps=T, cf ticket #1636) 196 !! but ensures reproductible results 197 !! with previous versions using split-explicit free surface 196 198 IF( ln_zps .AND. .NOT. ln_isfcav) & ! Partial steps: bottom before horizontal gradient 197 199 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! of t, s, rd at the last ocean level … … 201 203 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 202 204 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 203 204 ua(:,:,:) = 0._wp ! set dynamics trends to zero 205 va(:,:,:) = 0._wp 206 IF( lk_asminc .AND. ln_asmiau .AND. & 207 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 208 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 209 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 210 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 211 CALL dyn_ldf ( kstp ) ! lateral mixing 212 #if defined key_agrif 213 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 214 #endif 215 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 216 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 217 218 ua_sv(:,:,:) = ua(:,:,:) ! Save trends (barotropic trend has been fully updated at this stage) 219 va_sv(:,:,:) = va(:,:,:) 220 221 CALL div_hor( kstp ) ! Horizontal divergence (2nd call in time-split case) 222 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 223 CALL wzv ( kstp ) ! now cross-level velocity 224 ENDIF 225 226 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 227 ! diagnostics and outputs (ua, va, tsa used as workspace) 228 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 229 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 230 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 231 IF(.NOT.ln_cpl ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 232 IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports 233 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag 234 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 235 CALL dia_wri( kstp ) ! ocean model: outputs 236 ! 237 IF( ln_crs ) CALL crs_fld( kstp ) ! ocean model: online field coarsening & output 205 !!jc: fs simplification 206 207 ua(:,:,:) = 0._wp ! set dynamics trends to zero 208 va(:,:,:) = 0._wp 209 210 IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) & 211 CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 212 IF( lk_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends 213 #if defined key_agrif 214 IF(.NOT. Agrif_Root()) & 215 & CALL Agrif_Sponge_dyn ! momentum sponge 216 #endif 217 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 218 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 219 CALL dyn_ldf ( kstp ) ! lateral mixing 220 CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure 221 CALL dyn_spg ( kstp ) ! surface pressure gradient 222 223 ! With split-explicit free surface, since now transports have been updated and ssha as well 224 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 225 CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case) 226 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 227 CALL wzv ( kstp ) ! now cross-level velocity 228 ENDIF 229 230 CALL dyn_bfr ( kstp ) ! bottom friction 231 CALL dyn_zdf ( kstp ) ! vertical diffusion 232 233 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 234 ! diagnostics and outputs 235 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 236 IF( lk_floats ) CALL flo_stp ( kstp ) ! drifting Floats 237 IF( lk_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 238 IF(.NOT.ln_cpl ) CALL dia_fwb ( kstp ) ! Fresh water budget diagnostics 239 IF( lk_diadct ) CALL dia_dct ( kstp ) ! Transports 240 IF( lk_diaar5 ) CALL dia_ar5 ( kstp ) ! ar5 diag 241 IF( lk_diaharm ) CALL dia_harm ( kstp ) ! Tidal harmonic analysis 242 CALL dia_wri ( kstp ) ! ocean model: outputs 243 ! 244 IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output 238 245 239 246 #if defined key_top … … 241 248 ! Passive Tracer Model 242 249 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 243 CALL trc_stp ( kstp )! time-stepping244 #endif 245 246 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 247 ! Active tracers (ua, va used as workspace)248 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 249 tsa(:,:,:,:) = 0._wp! set tracer trends to zero250 CALL trc_stp ( kstp ) ! time-stepping 251 #endif 252 253 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 254 ! Active tracers 255 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 256 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 250 257 251 258 IF( lk_asminc .AND. ln_asmiau .AND. & 252 & ln_trainc ) CALL tra_asm_inc( kstp ) ! apply tracer assimilation increment 253 CALL tra_sbc ( kstp ) ! surface boundary condition 254 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 255 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 256 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 257 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 258 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 259 CALL tra_adv ( kstp ) ! horizontal & vertical advection 260 CALL tra_ldf ( kstp ) ! lateral mixing 259 & ln_trainc ) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment 260 CALL tra_sbc ( kstp ) ! surface boundary condition 261 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 262 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 263 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 264 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 265 IF( lk_bdy ) CALL bdy_tra_dmp ( kstp ) ! bdy damping trends 266 #if defined key_agrif 267 IF(.NOT. Agrif_Root()) & 268 & CALL Agrif_Sponge_tra ! tracers sponge 269 #endif 270 CALL tra_adv ( kstp ) ! horizontal & vertical advection 271 CALL tra_ldf ( kstp ) ! lateral mixing 261 272 262 273 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 263 IF( ln_diaptr ) CALL dia_ptr! Poleward adv/ldf TRansports diagnostics274 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 264 275 !!gm 265 266 #if defined key_agrif 267 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 268 #endif 269 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 270 271 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 272 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 273 CALL tra_nxt( kstp ) ! tracer fields at next time step 274 !!gm : why again a call to sto_pts ??? 275 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 276 !!gm 277 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 278 IF( ln_zps .AND. .NOT. ln_isfcav) & 279 & CALL zps_hde ( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps: before horizontal gradient 280 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 281 IF( ln_zps .AND. ln_isfcav) & 282 & CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top/bottom cells 283 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 284 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 285 ELSE ! centered hpg (eos then time stepping) 286 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 287 !!gm : why again a call to sto_pts ??? 288 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 289 !!gm 290 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 291 IF( ln_zps .AND. .NOT. ln_isfcav) & 292 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: bottom before horizontal gradient 293 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 294 IF( ln_zps .AND. ln_isfcav) & 295 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top/bottom cells 296 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 297 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 298 ENDIF 299 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 300 CALL tra_nxt( kstp ) ! tracer fields at next time step 301 ENDIF 302 303 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 304 ! Dynamics (tsa used as workspace) 305 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 306 IF( lk_dynspg_ts ) THEN 307 ! revert to previously computed momentum tendencies 308 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 309 ua(:,:,:) = ua_sv(:,:,:) 310 va(:,:,:) = va_sv(:,:,:) 311 312 CALL dyn_bfr( kstp ) ! bottom friction 313 CALL dyn_zdf( kstp ) ! vertical diffusion 314 ELSE 315 ua(:,:,:) = 0._wp ! set dynamics trends to zero 316 va(:,:,:) = 0._wp 317 318 IF( lk_asminc .AND. ln_asmiau .AND. & 319 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 320 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 321 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 322 CALL dyn_adv( kstp ) ! advection (vector or flux form) 323 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 324 CALL dyn_ldf( kstp ) ! lateral mixing 325 #if defined key_agrif 326 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 327 #endif 328 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 329 CALL dyn_bfr( kstp ) ! bottom friction 330 CALL dyn_zdf( kstp ) ! vertical diffusion 331 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 332 ENDIF 333 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 334 335 CALL ssh_swp( kstp ) ! swap of sea surface height 336 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 276 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 277 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 278 279 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 280 ! Set boundary conditions and Swap 281 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 282 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap 283 !! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields. 284 !! If so: 285 !! (i) no need to call agrif update at initialization time 286 !! (ii) no need to update "before" fields 287 !! 288 !! Apart from creating new tra_swp/dyn_swp routines, this however: 289 !! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between 290 !! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation", 291 !! e.g. a shift of the feedback interface inside child domain. 292 !! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 293 !! place. 294 !! 295 !!jc2: dynnxt must be the latest call. fse3t_b are indeed updated in that routine 296 CALL tra_nxt ( kstp ) ! finalize (bcs) tracer fields at next time step and swap 297 CALL dyn_nxt ( kstp ) ! finalize (bcs) velocities at next time step and swap 298 CALL ssh_swp ( kstp ) ! swap of sea surface height 299 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 337 300 ! 338 301 339 302 !!gm : This does not only concern the dynamics ==>>> add a new title 340 303 !!gm2: why ouput restart before AGRIF update? 341 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file 304 !! 305 !!jc: That would be better, but see comment above 306 !! 307 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 342 308 343 309 #if defined key_agrif … … 345 311 ! AGRIF 346 312 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 347 CALL Agrif_Integrate_ChildGrids( stp ) 348 349 IF ( Agrif_NbStepint().EQ.0 ) THEN 350 CALL Agrif_Update_Tra() ! Update active tracers 351 CALL Agrif_Update_Dyn() ! Update momentum 352 ENDIF 353 #endif 354 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 355 IF( lk_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 313 CALL Agrif_Integrate_ChildGrids( stp ) 314 315 IF ( Agrif_NbStepint().EQ.0 ) THEN ! AGRIF Update 316 !!jc in fact update i useless at last time step, but do it for global diagnostics 317 CALL Agrif_Update_Tra() ! Update active tracers 318 CALL Agrif_Update_Dyn() ! Update momentum 319 ENDIF 320 #endif 321 IF( ln_diahsb ) CALL dia_hsb ( kstp ) ! - ML - global conservation diagnostics 322 IF( lk_diaobs ) CALL dia_obs ( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 356 323 357 324 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 358 325 ! Control 359 326 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 360 CALL stp_ctl( kstp, indic )361 IF( indic < 0 362 363 364 ENDIF 365 IF( kstp == nit000 327 CALL stp_ctl ( kstp, indic ) 328 IF( indic < 0 ) THEN 329 CALL ctl_stop( 'step: indic < 0' ) 330 CALL dia_wri_state( 'output.abort', kstp ) 331 ENDIF 332 IF( kstp == nit000 ) THEN 366 333 CALL iom_close( numror ) ! close input ocean restart file 367 334 IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce … … 374 341 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 375 342 !!gm why lk_oasis and not lk_cpl ???? 376 IF( lk_oasis 343 IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 377 344 ! 378 345 #if defined key_iomput -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r5901 r6079 40 40 USE dynldf ! lateral momentum diffusion (dyn_ldf routine) 41 41 USE dynzdf ! vertical diffusion (dyn_zdf routine) 42 USE dynspg_oce ! surface pressure gradient (dyn_spg routine)43 42 USE dynspg ! surface pressure gradient (dyn_spg routine) 44 43 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
r3294 r6079 16 16 USE oce ! ocean dynamics and tracers variables 17 17 USE dom_oce ! ocean space and time domain variables 18 USE sol_oce ! ocean space and time domain variables19 18 USE in_out_manager ! I/O manager 20 19 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 21 20 USE lib_mpp ! distributed memory computing 22 USE dynspg_oce ! pressure gradient schemes23 21 USE c1d ! 1D vertical configuration 24 22 … … 43 41 !! ** Method : - Save the time step in numstp 44 42 !! - Print it each 50 time steps 45 !! - Print solver statistics in numsol 46 !! - Stop the run IF problem for the solver ( indec < 0 ) 43 !! - Stop the run IF problem ( indic < 0 ) 47 44 !! 48 45 !! ** Actions : 'time.step' file containing the last ocean time-step … … 50 47 !!---------------------------------------------------------------------- 51 48 INTEGER, INTENT( in ) :: kt ! ocean time-step index 52 INTEGER, INTENT( inout ) :: kindic ! indicator of solver convergence49 INTEGER, INTENT( inout ) :: kindic ! error indicator 53 50 !! 54 51 INTEGER :: ji, jj, jk ! dummy loop indices … … 143 140 IF( lk_c1d ) RETURN ! No log file in case of 1D vertical configuration 144 141 145 ! log file (solver or ssh statistics) 146 ! -------- 147 IF( lk_dynspg_flt ) THEN ! elliptic solver statistics (if required) 148 ! 149 IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps ! Solver 150 ! 151 IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN ! create a abort file if problem found 152 IF(lwp) THEN 153 WRITE(numout,*) ' stpctl: the elliptic solver DO not converge or explode' 154 WRITE(numout,*) ' ====== ' 155 WRITE(numout,9200) kt, niter, res, sqrt(epsr)/eps 156 WRITE(numout,*) 157 WRITE(numout,*) ' stpctl: output of last fields' 158 WRITE(numout,*) ' ====== ' 159 ENDIF 160 ENDIF 161 ! 162 ELSE !* ssh statistics (and others...) 163 IF( kt == nit000 .AND. lwp ) THEN ! open ssh statistics file (put in solver.stat file) 164 CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 165 ENDIF 166 ! 167 zssh2 = SUM( sshn(:,:) * sshn(:,:) * tmask_i(:,:) ) 168 IF( lk_mpp ) CALL mpp_sum( zssh2 ) ! sum over the global domain 169 ! 170 IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin ! ssh statistics 171 ! 142 ! log file (ssh statistics) 143 ! -------- !* ssh statistics (and others...) 144 IF( kt == nit000 .AND. lwp ) THEN ! open ssh statistics file (put in solver.stat file) 145 CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 172 146 ENDIF 147 ! 148 zssh2 = SUM( sshn(:,:) * sshn(:,:) * tmask_i(:,:) ) 149 IF( lk_mpp ) CALL mpp_sum( zssh2 ) ! sum over the global domain 150 ! 151 IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin ! ssh statistics 152 ! 173 153 174 154 9200 FORMAT('it:', i8, ' iter:', i4, ' r: ',e16.10, ' b: ',e16.10 ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/SAS_SRC
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/SAS_SRC/diawri.F90
r5901 r6079 26 26 USE dom_oce ! ocean space and time domain 27 27 USE zdf_oce ! ocean vertical physics 28 USE sol_oce ! solver variables29 28 USE sbc_oce ! Surface boundary condition: ocean fields 30 29 USE sbc_ice ! Surface boundary condition: ice fields -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/TOP_SRC
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/TOP_SRC/trcsub.F90
r5901 r6079 502 502 z1_rau0 = 0.5 / rau0 503 503 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 504 #if ! defined key_dynspg_ts 504 505 505 ! These lines are not necessary with time splitting since 506 506 ! boundary condition on sea level is set during ts loop … … 512 512 CALL lbc_lnk( ssha, 'T', 1. ) 513 513 #endif 514 #endif515 516 514 517 515 ! !------------------------------!
Note: See TracChangeset
for help on using the changeset viewer.