- Timestamp:
- 2017-06-01T13:24:36+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r6486 r8112 44 44 USE wrk_nemo ! Memory Allocation 45 45 USE timing ! Timing 46 USE iom 46 47 47 48 IMPLICIT NONE … … 379 380 !! 380 381 INTEGER :: ji, jj, jk ! dummy loop indices 381 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 382 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 383 !!---------------------------------------------------------------------- 384 ! 385 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 382 REAL(wp) :: zcoef0, zuap, zvap, znad, zvnum ! temporary scalars 383 REAL(wp), POINTER, DIMENSION(:,:) :: zunum, zudenom ! temporary numerators for u and v velocities (for debugging) 384 REAL(wp), POINTER, DIMENSION(:,:) :: zpb, zpt ! pressure at top and bottom of cells 385 REAL(wp), POINTER, DIMENSION(:,:) :: zzb, zzt ! geopotential height (gz) at top and bottom of cells 386 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj, zp3d, zz3d 387 !!---------------------------------------------------------------------- 388 ! 389 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zp3d, zz3d ) 390 CALL wrk_alloc( jpi,jpj, zunum, zudenom, zpb, zpt, zzb, zzt ) 386 391 ! 387 392 IF( kt == nit000 ) THEN … … 398 403 ENDIF 399 404 400 ! Surface value 401 DO jj = 2, jpjm1 402 DO ji = fs_2, fs_jpim1 ! vector opt. 403 ! hydrostatic pressure gradient along s-surfaces 404 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 405 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 406 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 407 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 408 ! s-coordinate pressure gradient correction 409 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 410 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 411 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 412 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 413 ! add to the general momentum trend 414 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 415 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 416 END DO 417 END DO 418 419 ! interior value (2=<jk=<jpkm1) 420 DO jk = 2, jpkm1 405 ! set the pressure values to zero at the top of the model (the atmospheric surface pressure is included in the barotropic step ?) 406 zpb(:,:) = 0._wp 407 zzb(:,:) = grav * sshn(:,:) ! surface height had been incorrectly set to zero; corrected on 24 Apr 2017 408 409 if(lwp) write(numout,*) 'hpg_sco: grav is ',grav 410 411 DO jk = 1, jpkm1 412 413 zpt(:,:) = zpb(:,:) 414 zzt(:,:) = zzb(:,:) 415 416 zp3d(:,:,jk) = zpt(:,:) 417 zz3d(:,:,jk) = zzt(:,:) 418 419 ! integrate the pressure from the top of the cell to the bottom of the cell; a factor rho_0 has been removed from the pressure 420 ! corrected to calculate for loops from 1 to jpj and 1 to jpi 25 Apr 2017 421 DO jj = 1, jpj 422 DO ji = 1, jpi ! vector opt. 423 !!$ zpb(ji,jj) = zpt(ji,jj) + grav * fse3w_n(ji,jj,jk) * ( znad + rhd(ji,jj,jk) ) 424 zpb(ji,jj) = zpt(ji,jj) + grav * fse3t_n(ji,jj,jk) * ( znad + rhd(ji,jj,jk) ) 425 zzb(ji,jj) = grav * ( sshn(ji,jj) - fsdepw_n(ji,jj,jk+1) ) ! surface height had been set to zero; corrected on 24 Apr 2017 426 END DO 427 END DO 428 429 ! calculate the Jacobian of the pressure forces 421 430 DO jj = 2, jpjm1 422 431 DO ji = fs_2, fs_jpim1 ! vector opt. 423 ! hydrostatic pressure gradient along s-surfaces 424 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 425 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 426 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 427 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 428 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 429 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 430 ! s-coordinate pressure gradient correction 431 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 432 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 433 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 434 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 435 ! add to the general momentum trend 436 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 437 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 432 zunum(ji, jj) = ( zpb(ji+1,jj) - zpt(ji , jj) )*( zzb(ji , jj) - zzt(ji+1, jj) ) - ( zpb(ji , jj) - zpt(ji+1, jj) )*( zzb(ji+1, jj) - zzt(ji , jj) ) 433 zudenom(ji, jj) = zpb(ji+1, jj) - zpt(ji , jj) + zpb(ji , jj) - zpt(ji+1, jj) 434 zhpi(ji, jj, jk) = zunum(ji, jj) * r1_e1u(ji, jj) / zudenom(ji, jj) 435 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 436 437 zvnum = ( zpb(ji,jj+1) - zpt(ji, jj ) )*( zzb(ji, jj ) - zzt(ji, jj+1) ) - ( zpb(ji, jj ) - zpt(ji, jj+1) )*( zzb(ji, jj+1) - zzt(ji, jj ) ) 438 zhpj(ji, jj, jk) = zvnum * r1_e2v(ji, jj) / ( zpb(ji, jj+1) - zpt(ji, jj ) + zpb(ji, jj ) - zpt(ji, jj+1) ) 439 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 438 440 END DO 439 441 END DO 440 END DO 441 ! 442 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 442 443 ! temporary diagnostic outputs (to be used on just a few timesteps on a small configuration) 444 !!$ DO jj = 1, jpj 445 !!$ write(numout, *) 'fse3w_n, jk, jj = ', jk, jj, (fse3w_n(ji,jj,jk), ji = 1, jpi) 446 !!$ write(numout, *) 'rhd, jk, jj = ', jk, jj, (rhd(ji,jj,jk), ji = 1, jpi) 447 !!$ write(numout, *) 'gdepw_n, jk+1, jj = ', jk+1, jj, (gdepw_n(ji,jj,jk+1), ji = 1, jpim1) 448 !!$ write(numout, *) 'r1_e1u, jk, jj = ', jk, jj, (r1_e1u(ji,jj), ji = 1, jpim1) 449 !!$ write(numout, *) 'zpt, jk, jj = ', jk, jj, (zpt(ji,jj), ji = 1, jpi) 450 !!$ write(numout, *) 'zpb, jk, jj = ', jk, jj, (zpb(ji,jj), ji = 1, jpi) 451 !!$ write(numout, *) 'zzt, jk, jj = ', jk, jj, (zzt(ji,jj), ji = 1, jpi) 452 !!$ write(numout, *) 'zzb, jk, jj = ', jk, jj, (zzb(ji,jj), ji = 1, jpi) 453 !!$ write(numout, *) 'zunum, jk, jj = ', jk, jj, (zunum(ji,jj), ji = 1, jpi) 454 !!$ write(numout, *) 'zudenom, jk, jj = ', jk, jj, (zudenom(ji,jj), ji = 1, jpi) 455 !!$ write(numout, *) 'tmask, jk, jj = ', jk, jj, (tmask(ji,jj,jk), ji = 1, jpi) 456 !!$ write(numout, *) 'umask, jk, jj = ', jk, jj, (umask(ji,jj,jk), ji = 1, jpi) 457 !!$ write(numout, *) 'vmask, jk, jj = ', jk, jj, (vmask(ji,jj,jk), ji = 1, jpi) 458 !!$ 459 !!$ write(numout, *) 'zhpi, jk, jj = ', jk, jj, (zhpi(ji,jj,jk), ji = 1, jpi) 460 !!$ write(numout, *) 'zhpj, jk, jj = ', jk, jj, (zhpj(ji,jj,jk), ji = 1, jpi) 461 !!$ END DO 462 463 END DO 464 ! 465 CALL iom_put("pressure" , zp3d(:,:,:) ) 466 CALL iom_put("hpgdepth" , zz3d(:,:,:) ) 467 CALL iom_put("rho" , rhd(:,:,:) ) 468 CALL iom_put("hpgi", zhpi(:,:,:) ) 469 CALL iom_put("hpgj", zhpj(:,:,:) ) 470 ! 471 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zp3d, zz3d ) 472 CALL wrk_dealloc( jpi,jpj, zpb, zpt, zzb, zzt ) 443 473 ! 444 474 END SUBROUTINE hpg_sco
Note: See TracChangeset
for help on using the changeset viewer.