Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r2715 r3294 14 14 !! - ! 2005-11 (G. Madec) style & small optimisation 15 15 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 !! ! (A. Coward) suppression of hel, wdj and rot options 16 18 !!---------------------------------------------------------------------- 17 19 … … 23 25 !! hpg_zps : z-coordinate plus partial steps (interpolation) 24 26 !! hpg_sco : s-coordinate (standard jacobian formulation) 25 !! hpg_hel : s-coordinate (helsinki modification)26 !! hpg_wdj : s-coordinate (weighted density jacobian)27 27 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 !! hpg_ rot : s-coordinate (ROTated axes scheme)28 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) 29 29 !!---------------------------------------------------------------------- 30 30 USE oce ! ocean dynamics and tracers … … 37 37 USE lbclnk ! lateral boundary condition 38 38 USE lib_mpp ! MPP library 39 USE wrk_nemo ! Memory Allocation 40 USE timing ! Timing 39 41 40 42 IMPLICIT NONE … … 48 50 LOGICAL , PUBLIC :: ln_hpg_zps = .FALSE. !: z-coordinate - partial steps (interpolation) 49 51 LOGICAL , PUBLIC :: ln_hpg_sco = .FALSE. !: s-coordinate (standard jacobian formulation) 50 LOGICAL , PUBLIC :: ln_hpg_hel = .FALSE. !: s-coordinate (helsinki modification)51 LOGICAL , PUBLIC :: ln_hpg_wdj = .FALSE. !: s-coordinate (weighted density jacobian)52 52 LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial) 53 LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme) 54 REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient 53 LOGICAL , PUBLIC :: ln_hpg_prj = .FALSE. !: s-coordinate (Pressure Jacobian scheme) 55 54 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 56 55 57 INTEGER :: nhpg = 0 ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)56 INTEGER :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 58 57 59 58 !! * Substitutions … … 77 76 !! - Save the trend (l_trddyn=T) 78 77 !!---------------------------------------------------------------------- 79 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released80 USE wrk_nemo, ONLY: ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2 ! 3D workspace81 !!82 78 INTEGER, INTENT(in) :: kt ! ocean time-step index 83 !!---------------------------------------------------------------------- 84 ! 85 IF( wrk_in_use(3, 1,2) ) THEN 86 CALL ctl_stop('dyn_hpg: requested workspace arrays are unavailable') ; RETURN 87 ENDIF 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 80 !!---------------------------------------------------------------------- 81 ! 82 IF( nn_timing == 1 ) CALL timing_start('dyn_hpg') 88 83 ! 89 84 IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn) 85 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 90 86 ztrdu(:,:,:) = ua(:,:,:) 91 87 ztrdv(:,:,:) = va(:,:,:) 92 88 ENDIF 93 89 ! 94 SELECT CASE ( nhpg ) ! Hydr astatic pressure gradient computation90 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation 95 91 CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate 96 92 CASE ( 1 ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation) 97 93 CASE ( 2 ) ; CALL hpg_sco ( kt ) ! s-coordinate (standard jacobian formulation) 98 CASE ( 3 ) ; CALL hpg_hel ( kt ) ! s-coordinate (helsinki modification) 99 CASE ( 4 ) ; CALL hpg_wdj ( kt ) ! s-coordinate (weighted density jacobian) 100 CASE ( 5 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 101 CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme) 94 CASE ( 3 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 95 CASE ( 4 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 102 96 END SELECT 103 97 ! … … 106 100 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 107 101 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt ) 102 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 108 103 ENDIF 109 104 ! … … 111 106 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 112 107 ! 113 IF( wrk_not_released(3, 1,2) ) CALL ctl_stop('dyn_hpg: failed to release workspace arrays')108 IF( nn_timing == 1 ) CALL timing_stop('dyn_hpg') 114 109 ! 115 110 END SUBROUTINE dyn_hpg … … 128 123 INTEGER :: ioptio = 0 ! temporary integer 129 124 !! 130 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel,&131 & ln_hpg_ wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma, ln_dynhpg_imp125 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 126 & ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp 132 127 !!---------------------------------------------------------------------- 133 128 ! … … 143 138 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 144 139 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 145 WRITE(numout,*) ' s-coord. (helsinki modification) ln_hpg_hel = ', ln_hpg_hel146 WRITE(numout,*) ' s-coord. (weighted density jacobian) ln_hpg_wdj = ', ln_hpg_wdj147 140 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 148 WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot 149 WRITE(numout,*) ' weighting coeff. (wdj scheme) rn_gamma = ', rn_gamma 141 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 150 142 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp 151 143 ENDIF 152 144 ! 153 IF( lk_vvl .AND. .NOT. ln_hpg_sco ) & 154 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco') 145 IF( ln_hpg_djc ) & 146 & CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method & 147 & currently disabled (bugs under investigation). Please select & 148 & either ln_hpg_sco or ln_hpg_prj instead') 149 ! 150 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) ) & 151 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 152 & the standard jacobian formulation hpg_sco or & 153 & the pressure jacobian formulation hpg_prj') 155 154 ! 156 155 ! ! Set nhpg from ln_hpg_... flags … … 158 157 IF( ln_hpg_zps ) nhpg = 1 159 158 IF( ln_hpg_sco ) nhpg = 2 160 IF( ln_hpg_hel ) nhpg = 3 161 IF( ln_hpg_wdj ) nhpg = 4 162 IF( ln_hpg_djc ) nhpg = 5 163 IF( ln_hpg_rot ) nhpg = 6 164 ! 165 ! ! Consitency check 159 IF( ln_hpg_djc ) nhpg = 3 160 IF( ln_hpg_prj ) nhpg = 4 161 ! 162 ! ! Consistency check 166 163 ioptio = 0 167 164 IF( ln_hpg_zco ) ioptio = ioptio + 1 168 165 IF( ln_hpg_zps ) ioptio = ioptio + 1 169 166 IF( ln_hpg_sco ) ioptio = ioptio + 1 170 IF( ln_hpg_hel ) ioptio = ioptio + 1171 IF( ln_hpg_wdj ) ioptio = ioptio + 1172 167 IF( ln_hpg_djc ) ioptio = ioptio + 1 173 IF( ln_hpg_ rot) ioptio = ioptio + 1168 IF( ln_hpg_prj ) ioptio = ioptio + 1 174 169 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 175 170 ! … … 193 188 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 194 189 !!---------------------------------------------------------------------- 195 USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace196 !!197 190 INTEGER, INTENT(in) :: kt ! ocean time-step index 198 191 !! 199 192 INTEGER :: ji, jj, jk ! dummy loop indices 200 193 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars 201 !!---------------------------------------------------------------------- 202 194 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 195 !!---------------------------------------------------------------------- 196 ! 197 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 198 ! 203 199 IF( kt == nit000 ) THEN 204 200 IF(lwp) WRITE(numout,*) … … 221 217 END DO 222 218 END DO 219 223 220 ! 224 221 ! interior value (2=<jk=<jpkm1) … … 242 239 END DO 243 240 ! 241 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 242 ! 244 243 END SUBROUTINE hpg_zco 245 244 … … 253 252 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 254 253 !!---------------------------------------------------------------------- 255 USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace256 !!257 254 INTEGER, INTENT(in) :: kt ! ocean time-step index 258 255 !! … … 260 257 INTEGER :: iku, ikv ! temporary integers 261 258 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 262 !!---------------------------------------------------------------------- 263 259 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 260 !!---------------------------------------------------------------------- 261 ! 262 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 263 ! 264 264 IF( kt == nit000 ) THEN 265 265 IF(lwp) WRITE(numout,*) … … 267 267 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization' 268 268 ENDIF 269 269 270 270 271 ! Local constant initialization … … 284 285 END DO 285 286 287 286 288 ! interior value (2=<jk=<jpkm1) 287 289 DO jk = 2, jpkm1 … … 303 305 END DO 304 306 END DO 307 305 308 306 309 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) … … 333 336 END DO 334 337 ! 338 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 339 ! 335 340 END SUBROUTINE hpg_zps 336 341 … … 354 359 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 355 360 !!---------------------------------------------------------------------- 356 USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace357 !!358 361 INTEGER, INTENT(in) :: kt ! ocean time-step index 359 362 !! 360 363 INTEGER :: ji, jj, jk ! dummy loop indices 361 364 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 362 !!---------------------------------------------------------------------- 363 365 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 366 !!---------------------------------------------------------------------- 367 ! 368 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 369 ! 364 370 IF( kt == nit000 ) THEN 365 371 IF(lwp) WRITE(numout,*) … … 417 423 END DO 418 424 ! 425 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 426 ! 419 427 END SUBROUTINE hpg_sco 420 421 422 SUBROUTINE hpg_hel( kt )423 !!---------------------------------------------------------------------424 !! *** ROUTINE hpg_hel ***425 !!426 !! ** Method : s-coordinate case.427 !! The now hydrostatic pressure gradient at a given level428 !! jk is computed by taking the vertical integral of the in-situ429 !! density gradient along the model level from the suface to that430 !! level. s-coordinates (ln_sco): a corrective term is added431 !! to the horizontal pressure gradient :432 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]433 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]434 !! add it to the general momentum trend (ua,va).435 !! ua = ua - 1/e1u * zhpi436 !! va = va - 1/e2v * zhpj437 !!438 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend439 !! - Save the trend (l_trddyn=T)440 !!----------------------------------------------------------------------441 USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace442 !!443 INTEGER, INTENT(in) :: kt ! ocean time-step index444 !!445 INTEGER :: ji, jj, jk ! dummy loop indices446 REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars447 !!----------------------------------------------------------------------448 449 IF( kt == nit000 ) THEN450 IF(lwp) WRITE(numout,*)451 IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend'452 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, helsinki modified scheme'453 ENDIF454 455 ! Local constant initialization456 zcoef0 = - grav * 0.5_wp457 458 ! Surface value459 DO jj = 2, jpjm1460 DO ji = fs_2, fs_jpim1 ! vector opt.461 ! hydrostatic pressure gradient along s-surfaces462 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) &463 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )464 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) &465 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )466 ! s-coordinate pressure gradient correction467 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &468 & * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)469 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &470 & * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)471 ! add to the general momentum trend472 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap473 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap474 END DO475 END DO476 !477 ! interior value (2=<jk=<jpkm1)478 DO jk = 2, jpkm1479 DO jj = 2, jpjm1480 DO ji = fs_2, fs_jpim1 ! vector opt.481 ! hydrostatic pressure gradient along s-surfaces482 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &483 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk) &484 & -fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk) ) &485 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) &486 & -fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) )487 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &488 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) &489 & -fse3t(ji,jj ,jk ) * rhd(ji,jj, jk) ) &490 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &491 & -fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) )492 ! s-coordinate pressure gradient correction493 zuap = - zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &494 & * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)495 zvap = - zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &496 & * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)497 ! add to the general momentum trend498 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap499 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap500 END DO501 END DO502 END DO503 !504 END SUBROUTINE hpg_hel505 506 507 SUBROUTINE hpg_wdj( kt )508 !!---------------------------------------------------------------------509 !! *** ROUTINE hpg_wdj ***510 !!511 !! ** Method : Weighted Density Jacobian (wdj) scheme (song 1998)512 !! The weighting coefficients from the namelist parameter rn_gamma513 !! (alpha=0.5-rn_gamma ; beta=1-alpha=0.5+rn_gamma514 !!515 !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998.516 !!----------------------------------------------------------------------517 USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace518 !!519 INTEGER, INTENT(in) :: kt ! ocean time-step index520 !!521 INTEGER :: ji, jj, jk ! dummy loop indices522 REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars523 REAL(wp) :: zalph , zbeta ! " "524 !!----------------------------------------------------------------------525 526 IF( kt == nit000 ) THEN527 IF(lwp) WRITE(numout,*)528 IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend'529 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ Weighted Density Jacobian'530 ENDIF531 532 ! Local constant initialization533 zcoef0 = - grav * 0.5_wp534 zalph = 0.5_wp - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma535 zbeta = 0.5_wp + rn_gamma ! (beta =1-alpha=0.5+rn_gamma536 537 ! Surface value (no ponderation)538 DO jj = 2, jpjm1539 DO ji = fs_2, fs_jpim1 ! vector opt.540 ! hydrostatic pressure gradient along s-surfaces541 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * rhd(ji+1,jj ,1) &542 & - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) )543 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * rhd(ji ,jj+1,1) &544 & - fse3w(ji ,jj ,1) * rhd(ji, jj ,1) )545 ! s-coordinate pressure gradient correction546 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &547 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)548 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &549 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)550 ! add to the general momentum trend551 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap552 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap553 END DO554 END DO555 556 ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta)557 DO jk = 2, jpkm1558 DO jj = 2, jpjm1559 DO ji = fs_2, fs_jpim1 ! vector opt.560 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &561 & * ( ( fsde3w(ji+1,jj,jk ) + fsde3w(ji,jj,jk ) &562 & - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &563 & * ( zalph * ( rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &564 & + zbeta * ( rhd (ji+1,jj,jk ) - rhd (ji,jj,jk ) ) ) &565 & - ( rhd (ji+1,jj,jk ) + rhd (ji,jj,jk ) &566 & - rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &567 & * ( zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &568 & + zbeta * ( fsde3w(ji+1,jj,jk ) - fsde3w(ji,jj,jk ) ) ) )569 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &570 & * ( ( fsde3w(ji,jj+1,jk ) + fsde3w(ji,jj,jk ) &571 & - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &572 & * ( zalph * ( rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &573 & + zbeta * ( rhd (ji,jj+1,jk ) - rhd (ji,jj,jk ) ) ) &574 & - ( rhd (ji,jj+1,jk ) + rhd (ji,jj,jk ) &575 & - rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &576 & * ( zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &577 & + zbeta * ( fsde3w(ji,jj+1,jk ) - fsde3w(ji,jj,jk ) ) ) )578 ! add to the general momentum trend579 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)580 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)581 END DO582 END DO583 END DO584 !585 END SUBROUTINE hpg_wdj586 587 428 588 429 SUBROUTINE hpg_djc( kt ) … … 594 435 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 595 436 !!---------------------------------------------------------------------- 596 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released597 USE oce , ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace598 USE wrk_nemo, ONLY: drhox => wrk_3d_1 , dzx => wrk_3d_2599 USE wrk_nemo, ONLY: drhou => wrk_3d_3 , dzu => wrk_3d_4 , rho_i => wrk_3d_5600 USE wrk_nemo, ONLY: drhoy => wrk_3d_6 , dzy => wrk_3d_7601 USE wrk_nemo, ONLY: drhov => wrk_3d_8 , dzv => wrk_3d_9 , rho_j => wrk_3d_10602 USE wrk_nemo, ONLY: drhoz => wrk_3d_11 , dzz => wrk_3d_12603 USE wrk_nemo, ONLY: drhow => wrk_3d_13 , dzw => wrk_3d_14604 USE wrk_nemo, ONLY: rho_k => wrk_3d_15605 !!606 437 INTEGER, INTENT(in) :: kt ! ocean time-step index 607 438 !! … … 610 441 REAL(wp) :: z1_10, cffu, cffx ! " " 611 442 REAL(wp) :: z1_12, cffv, cffy ! " " 612 !!---------------------------------------------------------------------- 613 614 IF( wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) THEN 615 CALL ctl_stop('dyn:hpg_djc: requested workspace arrays unavailable') ; RETURN 616 ENDIF 443 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 444 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 445 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 446 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 447 !!---------------------------------------------------------------------- 448 ! 449 CALL wrk_alloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 450 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 451 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 452 ! 617 453 618 454 IF( kt == nit000 ) THEN … … 811 647 END DO 812 648 ! 813 IF( wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) & 814 CALL ctl_stop('dyn:hpg_djc: failed to release workspace arrays') 649 CALL wrk_dealloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 650 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 651 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 815 652 ! 816 653 END SUBROUTINE hpg_djc 817 654 818 655 819 SUBROUTINE hpg_ rot( kt )656 SUBROUTINE hpg_prj( kt ) 820 657 !!--------------------------------------------------------------------- 821 !! *** ROUTINE hpg_rot *** 822 !! 823 !! ** Method : rotated axes scheme (Thiem and Berntsen 2005) 824 !! 825 !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005. 826 !!---------------------------------------------------------------------- 827 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 828 USE oce , ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace 829 USE wrk_nemo, ONLY: zdistr => wrk_2d_1 , zsina => wrk_2d_2 , zcosa => wrk_2d_3 830 USE wrk_nemo, ONLY: zhpiorg => wrk_3d_1 , zhpirot => wrk_3d_2 831 USE wrk_nemo, ONLY: zhpitra => wrk_3d_3 , zhpine => wrk_3d_4 832 USE wrk_nemo, ONLY: zhpjorg => wrk_3d_5 , zhpjrot => wrk_3d_6 833 USE wrk_nemo, ONLY: zhpjtra => wrk_3d_7 , zhpjne => wrk_3d_8 834 !! 835 INTEGER, INTENT(in) :: kt ! ocean time-step index 836 !! 837 INTEGER :: ji, jj, jk ! dummy loop indices 838 REAL(wp) :: zforg, zcoef0, zuap, zmskd1, zmskd1m ! temporary scalar 839 REAL(wp) :: zfrot , zvap, zmskd2, zmskd2m ! " " 840 !!---------------------------------------------------------------------- 841 842 IF( wrk_in_use(2, 1,2,3) .OR. & 843 wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN 844 CALL ctl_stop('dyn:hpg_rot: requested workspace arrays unavailable') ; RETURN 845 ENDIF 846 658 !! *** ROUTINE hpg_prj *** 659 !! 660 !! ** Method : s-coordinate case. 661 !! A Pressure-Jacobian horizontal pressure gradient method 662 !! based on the constrained cubic-spline interpolation for 663 !! all vertical coordinate systems 664 !! 665 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 666 !! - Save the trend (l_trddyn=T) 667 !! 668 !!---------------------------------------------------------------------- 669 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 670 INTEGER, INTENT(in) :: kt ! ocean time-step index 671 !! 672 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 673 REAL(wp) :: zcoef0, znad ! temporary scalars 674 !! 675 !! The local variables for the correction term 676 INTEGER :: jk1, jis, jid, jjs, jjd 677 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 678 REAL(wp) :: zrhdt1 679 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 680 INTEGER :: zbhitwe, zbhitns 681 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdeptht, zrhh 682 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 683 !!---------------------------------------------------------------------- 684 ! 685 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 686 CALL wrk_alloc( jpi,jpj,jpk, zdeptht, zrhh ) 687 ! 847 688 IF( kt == nit000 ) THEN 848 689 IF(lwp) WRITE(numout,*) 849 IF(lwp) WRITE(numout,*) 'dyn:hpg_ rot: hydrostatic pressure gradient trend'850 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, rotated axes scheme used'690 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 691 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 851 692 ENDIF 852 693 853 ! ------------------------------- 854 ! Local constant initialization 855 ! ------------------------------- 856 zcoef0 = - grav * 0.5_wp 857 zforg = 0.95_wp 858 zfrot = 1._wp - zforg 859 860 ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance) 861 zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)*e1f(:,:) + e2f(:,:)*e1f(:,:) ) 862 863 ! sinus and cosinus of diagonal angle at F-point 864 zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) ) 865 zcosa(:,:) = COS( zsina(:,:) ) 866 zsina(:,:) = SIN( zsina(:,:) ) 867 868 ! --------------- 869 ! Surface value 870 ! --------------- 871 ! compute and add to the general trend the pressure gradients along the axes 872 DO jj = 2, jpjm1 873 DO ji = fs_2, fs_jpim1 ! vector opt. 874 ! hydrostatic pressure gradient along s-surfaces 875 zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,1) * rhd(ji+1,jj,1) & 876 & - fse3t(ji ,jj,1) * rhd(ji ,jj,1) ) 877 zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,1) * rhd(ji,jj+1,1) & 878 & - fse3t(ji,jj ,1) * rhd(ji,jj ,1) ) 879 ! s-coordinate pressure gradient correction 880 zuap = -zcoef0 * ( rhd (ji+1,jj ,1) + rhd (ji,jj,1) ) & 881 & * ( fsdept(ji+1,jj ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 882 zvap = -zcoef0 * ( rhd (ji ,jj+1,1) + rhd (ji,jj,1) ) & 883 & * ( fsdept(ji ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 884 ! add to the general momentum trend 885 ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap ) 886 va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap ) 887 END DO 888 END DO 889 890 ! compute the pressure gradients in the diagonal directions 891 DO jj = 1, jpjm1 892 DO ji = 1, fs_jpim1 ! vector opt. 893 zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji ,jj,1) ! mask in the 1st diagnonal 894 zmskd2 = tmask(ji ,jj+1,1) * tmask(ji+1,jj,1) ! mask in the 2nd diagnonal 895 ! hydrostatic pressure gradient along s-surfaces 896 zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1) & 897 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) ) 898 zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) & 899 & - fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) ) 900 ! s-coordinate pressure gradient correction 901 zuap = -zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,1) + rhd (ji ,jj,1) ) & 902 & * ( fsdept(ji+1,jj+1,1) - fsdept(ji ,jj,1) ) 903 zvap = -zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,1) + rhd (ji+1,jj,1) ) & 904 & * ( fsdept(ji ,jj+1,1) - fsdept(ji+1,jj,1) ) 905 ! back rotation 906 zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) & 907 & - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 908 zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) & 909 & + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 910 END DO 911 END DO 912 913 ! interpolate and add to the general trend the diagonal gradient 914 DO jj = 2, jpjm1 915 DO ji = fs_2, fs_jpim1 ! vector opt. 916 ! averaging 917 zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji ,jj-1,1) ) 918 zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj ,1) ) 919 ! add to the general momentum trend 920 ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1) 921 va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1) 922 END DO 923 END DO 924 925 ! ----------------- 926 ! 2. interior value (2=<jk=<jpkm1) 927 ! ----------------- 928 ! compute and add to the general trend the pressure gradients along the axes 929 DO jk = 2, jpkm1 930 DO jj = 2, jpjm1 931 DO ji = fs_2, fs_jpim1 ! vector opt. 932 ! hydrostatic pressure gradient along s-surfaces 933 zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1) & 934 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk ) & 935 & - fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk ) & 936 & + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) & 937 & - fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) ) 938 zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1) & 939 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk ) & 940 & - fse3t(ji,jj ,jk ) * rhd(ji,jj, jk ) & 941 & + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) & 942 & - fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) ) 943 ! s-coordinate pressure gradient correction 944 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) & 945 & * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 946 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) & 947 & * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 948 ! add to the general momentum trend 949 ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap ) 950 va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap ) 694 !!---------------------------------------------------------------------- 695 ! Local constant initialization 696 zcoef0 = - grav 697 znad = 0.0_wp 698 IF( lk_vvl ) znad = 1._wp 699 700 ! Clean 3-D work arrays 701 zhpi(:,:,:) = 0._wp 702 zrhh(:,:,:) = rhd(:,:,:) 703 704 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 705 DO jj = 1, jpj 706 DO ji = 1, jpi 707 jk = mbathy(ji,jj) 708 IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp 709 ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 710 ELSE IF(jk < jpkm1) THEN 711 DO jkk = jk+1, jpk 712 zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1), & 713 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 714 END DO 715 ENDIF 716 END DO 717 END DO 718 719 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)" 720 DO jj = 1, jpj 721 DO ji = 1, jpi 722 zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 723 zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad 724 DO jk = 2, jpk 725 zdeptht(ji,jj,jk) = zdeptht(ji,jj,jk-1) + fse3w(ji,jj,jk) 726 END DO 727 END DO 728 END DO 729 730 DO jk = 1, jpkm1 731 DO jj = 1, jpj 732 DO ji = 1, jpi 733 fsp(ji,jj,jk) = zrhh(ji,jj,jk) 734 xsp(ji,jj,jk) = zdeptht(ji,jj,jk) 735 END DO 736 END DO 737 END DO 738 739 ! Construct the vertical density profile with the 740 ! constrained cubic spline interpolation 741 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 742 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 743 744 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 745 DO jj = 2, jpj 746 DO ji = 2, jpi 747 zrhdt1 = zrhh(ji,jj,1) - interp3(zdeptht(ji,jj,1),asp(ji,jj,1), & 748 bsp(ji,jj,1), csp(ji,jj,1), & 749 dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1) 750 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 751 752 ! assuming linear profile across the top half surface layer 753 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1 754 END DO 755 END DO 756 757 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 758 DO jk = 2, jpkm1 759 DO jj = 2, jpj 760 DO ji = 2, jpi 761 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 762 integ2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),& 763 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), & 764 csp(ji,jj,jk-1), dsp(ji,jj,jk-1)) 765 END DO 766 END DO 767 END DO 768 769 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 770 DO jj = 2, jpjm1 771 DO ji = 2, jpim1 772 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 773 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 774 END DO 775 END DO 776 777 DO jk = 2, jpkm1 778 DO jj = 2, jpjm1 779 DO ji = 2, jpim1 780 zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 781 zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 782 END DO 783 END DO 784 END DO 785 786 DO jk = 1, jpkm1 787 DO jj = 2, jpjm1 788 DO ji = 2, jpim1 789 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 790 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 791 END DO 792 END DO 793 END DO 794 795 DO jk = 1, jpkm1 796 DO jj = 2, jpjm1 797 DO ji = 2, jpim1 798 zpwes = 0._wp; zpwed = 0._wp 799 zpnss = 0._wp; zpnsd = 0._wp 800 zuijk = zu(ji,jj,jk) 801 zvijk = zv(ji,jj,jk) 802 803 !!!!! for u equation 804 IF( jk <= mbku(ji,jj) ) THEN 805 IF( -zdeptht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN 806 jis = ji + 1; jid = ji 807 ELSE 808 jis = ji; jid = ji +1 809 ENDIF 810 811 ! integrate the pressure on the shallow side 812 jk1 = jk 813 zbhitwe = 0 814 DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk ) 815 IF( jk1 == mbku(ji,jj) ) THEN 816 zbhitwe = 1 817 EXIT 818 ENDIF 819 zdeps = MIN(zdeptht(jis,jj,jk1+1), -zuijk) 820 zpwes = zpwes + & 821 integ2(zdeptht(jis,jj,jk1), zdeps, & 822 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 823 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 824 jk1 = jk1 + 1 825 END DO 826 827 IF(zbhitwe == 1) THEN 828 zuijk = -zdeptht(jis,jj,jk1) 829 ENDIF 830 831 ! integrate the pressure on the deep side 832 jk1 = jk 833 zbhitwe = 0 834 DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk ) 835 IF( jk1 == 1 ) THEN 836 zbhitwe = 1 837 EXIT 838 ENDIF 839 zdeps = MAX(zdeptht(jid,jj,jk1-1), -zuijk) 840 zpwed = zpwed + & 841 integ2(zdeps, zdeptht(jid,jj,jk1), & 842 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 843 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 844 jk1 = jk1 - 1 845 END DO 846 847 IF( zbhitwe == 1 ) THEN 848 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 849 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), & 850 bsp(jid,jj,1), csp(jid,jj,1), & 851 dsp(jid,jj,1)) * zdeps 852 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 853 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 854 ENDIF 855 856 ! update the momentum trends in u direction 857 858 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 859 IF( lk_vvl ) THEN 860 zdpdx2 = zcoef0 / e1u(ji,jj) * & 861 ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 862 ELSE 863 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 864 ENDIF 865 866 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 867 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 868 ENDIF 869 870 !!!!! for v equation 871 IF( jk <= mbkv(ji,jj) ) THEN 872 IF( -zdeptht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN 873 jjs = jj + 1; jjd = jj 874 ELSE 875 jjs = jj ; jjd = jj + 1 876 ENDIF 877 878 ! integrate the pressure on the shallow side 879 jk1 = jk 880 zbhitns = 0 881 DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk ) 882 IF( jk1 == mbkv(ji,jj) ) THEN 883 zbhitns = 1 884 EXIT 885 ENDIF 886 zdeps = MIN(zdeptht(ji,jjs,jk1+1), -zvijk) 887 zpnss = zpnss + & 888 integ2(zdeptht(ji,jjs,jk1), zdeps, & 889 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 890 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 891 jk1 = jk1 + 1 892 END DO 893 894 IF(zbhitns == 1) THEN 895 zvijk = -zdeptht(ji,jjs,jk1) 896 ENDIF 897 898 ! integrate the pressure on the deep side 899 jk1 = jk 900 zbhitns = 0 901 DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk ) 902 IF( jk1 == 1 ) THEN 903 zbhitns = 1 904 EXIT 905 ENDIF 906 zdeps = MAX(zdeptht(ji,jjd,jk1-1), -zvijk) 907 zpnsd = zpnsd + & 908 integ2(zdeps, zdeptht(ji,jjd,jk1), & 909 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 910 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 911 jk1 = jk1 - 1 912 END DO 913 914 IF( zbhitns == 1 ) THEN 915 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 916 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), & 917 bsp(ji,jjd,1), csp(ji,jjd,1), & 918 dsp(ji,jjd,1) ) * zdeps 919 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 920 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 921 ENDIF 922 923 ! update the momentum trends in v direction 924 925 zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 926 IF( lk_vvl ) THEN 927 zdpdy2 = zcoef0 / e2v(ji,jj) * & 928 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 929 ELSE 930 zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 931 ENDIF 932 933 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 934 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 935 ENDIF 936 937 938 END DO 939 END DO 940 END DO 941 ! 942 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 943 CALL wrk_dealloc( jpi,jpj,jpk, zdeptht, zrhh ) 944 ! 945 END SUBROUTINE hpg_prj 946 947 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 948 !!---------------------------------------------------------------------- 949 !! *** ROUTINE cspline *** 950 !! 951 !! ** Purpose : constrained cubic spline interpolation 952 !! 953 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 954 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 955 !! 956 !!---------------------------------------------------------------------- 957 IMPLICIT NONE 958 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 959 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 960 ! the interpoated function 961 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 962 ! 2: Linear 963 964 ! Local Variables 965 INTEGER :: ji, jj, jk ! dummy loop indices 966 INTEGER :: jpi, jpj, jpkm1 967 REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 968 REAL(wp) :: zdxtmp1, zdxtmp2, zalpha 969 REAL(wp) :: zdf(size(fsp,3)) 970 !!---------------------------------------------------------------------- 971 972 jpi = size(fsp,1) 973 jpj = size(fsp,2) 974 jpkm1 = size(fsp,3) - 1 975 976 977 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 978 DO ji = 1, jpi 979 DO jj = 1, jpj 980 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 981 ! DO jk = 2, jpkm1-1 982 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 983 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 984 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 985 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 986 ! 987 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 988 ! 989 ! IF(zdf1 * zdf2 <= 0._wp) THEN 990 ! zdf(jk) = 0._wp 991 ! ELSE 992 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 993 ! ENDIF 994 ! END DO 995 996 !!Simply geometric average 997 DO jk = 2, jpkm1-1 998 zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1)) 999 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk)) 1000 1001 IF(zdf1 * zdf2 <= 0._wp) THEN 1002 zdf(jk) = 0._wp 1003 ELSE 1004 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1005 ENDIF 1006 END DO 1007 1008 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1009 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1010 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1011 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1012 & 0.5_wp * zdf(jpkm1 - 1) 1013 1014 DO jk = 1, jpkm1 - 1 1015 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1016 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1017 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1018 zddf1 = -2._wp * ztmp1 + ztmp2 1019 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1020 zddf2 = 2._wp * ztmp1 - ztmp2 1021 1022 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1023 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1024 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1025 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1026 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1027 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1028 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1029 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1030 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1031 END DO 951 1032 END DO 952 1033 END DO 953 END DO 954 955 ! compute the pressure gradients in the diagonal directions 956 DO jk = 2, jpkm1 957 DO jj = 1, jpjm1 958 DO ji = 1, fs_jpim1 ! vector opt. 959 zmskd1 = tmask(ji+1,jj+1,jk ) * tmask(ji ,jj,jk ) ! level jk mask in the 1st diagnonal 960 zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji ,jj,jk-1) ! level jk-1 " " 961 zmskd2 = tmask(ji ,jj+1,jk ) * tmask(ji+1,jj,jk ) ! level jk mask in the 2nd diagnonal 962 zmskd2m = tmask(ji ,jj+1,jk-1) * tmask(ji+1,jj,jk-1) ! level jk-1 " " 963 ! hydrostatic pressure gradient along s-surfaces 964 zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1) & 965 & + zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,jk ) * rhd(ji+1,jj+1,jk) & 966 & -fse3t(ji ,jj ,jk ) * rhd(ji ,jj ,jk) ) & 967 & + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1) & 968 & -fse3t(ji ,jj ,jk-1) * rhd(ji ,jj ,jk-1) ) 969 zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1) & 970 & + zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,jk ) * rhd(ji ,jj+1,jk) & 971 & -fse3t(ji+1,jj ,jk ) * rhd(ji+1,jj, jk) ) & 972 & + zdistr(ji,jj) * zmskd2m * ( fse3t(ji ,jj+1,jk-1) * rhd(ji ,jj+1,jk-1) & 973 & -fse3t(ji+1,jj ,jk-1) * rhd(ji+1,jj, jk-1) ) 974 ! s-coordinate pressure gradient correction 975 zuap = - zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,jk) + rhd (ji ,jj,jk) ) & 976 & * ( fsdept(ji+1,jj+1,jk) - fsdept(ji ,jj,jk) ) 977 zvap = - zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,jk) + rhd (ji+1,jj,jk) ) & 978 & * ( fsdept(ji ,jj+1,jk) - fsdept(ji+1,jj,jk) ) 979 ! back rotation 980 zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) & 981 & - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 982 zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) & 983 & + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 1034 1035 ELSE IF (polynomial_type == 2) THEN ! Linear 1036 DO ji = 1, jpi 1037 DO jj = 1, jpj 1038 DO jk = 1, jpkm1-1 1039 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1040 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1041 1042 dsp(ji,jj,jk) = 0._wp 1043 csp(ji,jj,jk) = 0._wp 1044 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1045 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1046 END DO 984 1047 END DO 985 1048 END DO 986 END DO 987 988 ! interpolate and add to the general trend 989 DO jk = 2, jpkm1 990 DO jj = 2, jpjm1 991 DO ji = fs_2, fs_jpim1 ! vector opt. 992 ! averaging 993 zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji ,jj-1,jk) ) 994 zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj ,jk) ) 995 ! add to the general momentum trend 996 ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk) 997 va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk) 998 END DO 999 END DO 1000 END DO 1001 ! 1002 IF( wrk_not_released(2, 1,2,3) .OR. & 1003 wrk_not_released(3, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:hpg_rot: failed to release workspace arrays') 1004 ! 1005 END SUBROUTINE hpg_rot 1049 1050 ELSE 1051 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1052 ENDIF 1053 1054 1055 END SUBROUTINE cspline 1056 1057 1058 FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f) 1059 !!---------------------------------------------------------------------- 1060 !! *** ROUTINE interp1 *** 1061 !! 1062 !! ** Purpose : 1-d linear interpolation 1063 !! 1064 !! ** Method : 1065 !! interpolation is straight forward 1066 !! extrapolation is also permitted (no value limit) 1067 !! 1068 !!---------------------------------------------------------------------- 1069 IMPLICIT NONE 1070 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1071 REAL(wp) :: f ! result of the interpolation (extrapolation) 1072 REAL(wp) :: zdeltx 1073 !!---------------------------------------------------------------------- 1074 1075 zdeltx = xr - xl 1076 IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 1077 f = 0.5_wp * (fl + fr) 1078 ELSE 1079 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1080 ENDIF 1081 1082 END FUNCTION interp1 1083 1084 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1085 !!---------------------------------------------------------------------- 1086 !! *** ROUTINE interp1 *** 1087 !! 1088 !! ** Purpose : 1-d constrained cubic spline interpolation 1089 !! 1090 !! ** Method : cubic spline interpolation 1091 !! 1092 !!---------------------------------------------------------------------- 1093 IMPLICIT NONE 1094 REAL(wp), INTENT(in) :: x, a, b, c, d 1095 REAL(wp) :: f ! value from the interpolation 1096 !!---------------------------------------------------------------------- 1097 1098 f = a + x* ( b + x * ( c + d * x ) ) 1099 1100 END FUNCTION interp2 1101 1102 1103 FUNCTION interp3(x, a, b, c, d) RESULT(f) 1104 !!---------------------------------------------------------------------- 1105 !! *** ROUTINE interp1 *** 1106 !! 1107 !! ** Purpose : Calculate the first order of deriavtive of 1108 !! a cubic spline function y=a+b*x+c*x^2+d*x^3 1109 !! 1110 !! ** Method : f=dy/dx=b+2*c*x+3*d*x^2 1111 !! 1112 !!---------------------------------------------------------------------- 1113 IMPLICIT NONE 1114 REAL(wp), INTENT(in) :: x, a, b, c, d 1115 REAL(wp) :: f ! value from the interpolation 1116 !!---------------------------------------------------------------------- 1117 1118 f = b + x * ( 2._wp * c + 3._wp * d * x) 1119 1120 END FUNCTION interp3 1121 1122 1123 FUNCTION integ2(xl, xr, a, b, c, d) RESULT(f) 1124 !!---------------------------------------------------------------------- 1125 !! *** ROUTINE interp1 *** 1126 !! 1127 !! ** Purpose : 1-d constrained cubic spline integration 1128 !! 1129 !! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr 1130 !! 1131 !!---------------------------------------------------------------------- 1132 IMPLICIT NONE 1133 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1134 REAL(wp) :: za1, za2, za3 1135 REAL(wp) :: f ! integration result 1136 !!---------------------------------------------------------------------- 1137 1138 za1 = 0.5_wp * b 1139 za2 = c / 3.0_wp 1140 za3 = 0.25_wp * d 1141 1142 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1143 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1144 1145 END FUNCTION integ2 1146 1006 1147 1007 1148 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
