Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Property svn:eol-style deleted
r2466 r2528 4 4 !! Ocean dynamics: hydrostatic pressure gradient trend 5 5 !!====================================================================== 6 !! History : 1.0 ! 87-09 (P. Andrich, M.-A. Foujols) hpg_zco: Original code 7 !! 5.0 ! 91-11 (G. Madec) 8 !! 7.0 ! 96-01 (G. Madec) hpg_sco: Original code for s-coordinates 9 !! 8.0 ! 97-05 (G. Madec) split dynber into dynkeg and dynhpg 10 !! 8.5 ! 02-07 (G. Madec) F90: Free form and module 11 !! 8.5 ! 02-08 (A. Bozec) hpg_zps: Original code 12 !! 9.0 ! 05-10 (A. Beckmann, B.W. An) various s-coordinate options 13 !! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot 14 !! 9.0 ! 05-11 (G. Madec) style & small optimisation 6 !! History : OPA ! 1987-09 (P. Andrich, M.-A. Foujols) hpg_zco: Original code 7 !! 5.0 ! 1991-11 (G. Madec) 8 !! 7.0 ! 1996-01 (G. Madec) hpg_sco: Original code for s-coordinates 9 !! 8.0 ! 1997-05 (G. Madec) split dynber into dynkeg and dynhpg 10 !! 8.5 ! 2002-07 (G. Madec) F90: Free form and module 11 !! 8.5 ! 2002-08 (A. Bozec) hpg_zps: Original code 12 !! NEMO 1.0 ! 2005-10 (A. Beckmann, B.W. An) various s-coordinate options 13 !! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot 14 !! - ! 2005-11 (G. Madec) style & small optimisation 15 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 15 16 !!---------------------------------------------------------------------- 16 17 … … 18 19 !! dyn_hpg : update the momentum trend with the now horizontal 19 20 !! gradient of the hydrostatic pressure 20 !! hpg_ctl: initialisation and control of options21 !! dyn_hpg_init : initialisation and control of options 21 22 !! hpg_zco : z-coordinate scheme 22 23 !! hpg_zps : z-coordinate plus partial steps (interpolation) … … 39 40 PRIVATE 40 41 41 PUBLIC dyn_hpg ! routine called by step module 42 PUBLIC dyn_hpg ! routine called by step module 43 PUBLIC dyn_hpg_init ! routine called by opa module 42 44 43 45 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient … … 49 51 LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial) 50 52 LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme) 51 REAL(wp), PUBLIC :: rn_gamma = 0. e0!: weighting coefficient53 REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient 52 54 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 53 55 … … 58 60 # include "vectopt_loop_substitute.h90" 59 61 !!---------------------------------------------------------------------- 60 !! OPA 9.0 , LOCEAN-IPSL (2005)61 !! $Id$ 62 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)62 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 63 !! $Id$ 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 63 65 !!---------------------------------------------------------------------- 64 65 66 CONTAINS 66 67 … … 79 80 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrdu, ztrdv ! 3D temporary workspace 80 81 !!---------------------------------------------------------------------- 81 82 IF( kt == nit000 ) CALL hpg_ctl ! initialisation & control of options 83 82 ! 84 83 IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn) 85 84 ztrdu(:,:,:) = ua(:,:,:) 86 85 ztrdv(:,:,:) = va(:,:,:) 87 86 ENDIF 88 87 ! 89 88 SELECT CASE ( nhpg ) ! Hydrastatic pressure gradient computation 90 89 CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate … … 96 95 CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme) 97 96 END SELECT 98 97 ! 99 98 IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 100 99 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) … … 109 108 110 109 111 SUBROUTINE hpg_ctl112 !!---------------------------------------------------------------------- 113 !! *** ROUTINE hpg_ctl***110 SUBROUTINE dyn_hpg_init 111 !!---------------------------------------------------------------------- 112 !! *** ROUTINE dyn_hpg_init *** 114 113 !! 115 114 !! ** Purpose : initializations for the hydrostatic pressure gradient … … 121 120 INTEGER :: ioptio = 0 ! temporary integer 122 121 !! 123 ! NAMELIST/namdyn_hpg/ ln_hpg_zco , ln_hpg_zps , ln_hpg_sco, ln_hpg_hel, & 124 ! & ln_hpg_wdj , ln_hpg_djc , ln_hpg_rot, rn_gamma , & 125 ! & ln_dynhpg_imp 126 !!---------------------------------------------------------------------- 127 128 ! REWIND ( numnam ) ! Namelist namdyn_hpg : already read in opa.F90 module 129 ! READ ( numnam, namdyn_hpg ) 130 131 IF(lwp) THEN ! Control print 122 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel, & 123 & ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma , ln_dynhpg_imp 124 !!---------------------------------------------------------------------- 125 ! 126 REWIND( numnam ) ! Read Namelist namdyn_hpg 127 READ ( numnam, namdyn_hpg ) 128 ! 129 IF(lwp) THEN ! Control print 132 130 WRITE(numout,*) 133 WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient'134 WRITE(numout,*) '~~~~~~~ '131 WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation' 132 WRITE(numout,*) '~~~~~~~~~~~~' 135 133 WRITE(numout,*) ' Namelist namdyn_hpg : choice of hpg scheme' 136 134 WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco … … 144 142 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp 145 143 ENDIF 146 147 IF( lk_vvl .AND. .NOT. ln_hpg_sco ) THEN 148 CALL ctl_stop( 'hpg_ctl : variable volume key_vvl compatible only with the standard jacobian formulation hpg_sco') 149 ENDIF 150 144 ! 145 IF( lk_vvl .AND. .NOT. ln_hpg_sco ) & 146 & CALL ctl_stop( 'dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco') 147 ! 151 148 ! ! Set nhpg from ln_hpg_... flags 152 149 IF( ln_hpg_zco ) nhpg = 0 … … 157 154 IF( ln_hpg_djc ) nhpg = 5 158 155 IF( ln_hpg_rot ) nhpg = 6 159 156 ! 160 157 ! ! Consitency check 161 158 ioptio = 0 … … 168 165 IF( ln_hpg_rot ) ioptio = ioptio + 1 169 166 IF ( ioptio /= 1 ) CALL ctl_stop( ' NO or several hydrostatic pressure gradient options used' ) 170 171 ! 172 END SUBROUTINE hpg_ctl 167 ! 168 END SUBROUTINE dyn_hpg_init 173 169 174 170 … … 205 201 206 202 ! Local constant initialization 207 zcoef0 = - grav * 0.5 203 zcoef0 = - grav * 0.5_wp 208 204 209 205 ! Surface value … … 268 264 269 265 ! Local constant initialization 270 zcoef0 = - grav * 0.5 271 272 ! Surface value 266 zcoef0 = - grav * 0.5_wp 267 268 ! Surface value (also valid in partial step case) 273 269 DO jj = 2, jpjm1 274 270 DO ji = fs_2, fs_jpim1 ! vector opt. … … 303 299 END DO 304 300 305 ! partial steps correction at the last level ( new gradient with intgrd.F)301 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) 306 302 # if defined key_vectopt_loop 307 303 jj = 1 … … 311 307 DO ji = 2, jpim1 312 308 # endif 313 iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1314 ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1309 iku = mbku(ji,jj) 310 ikv = mbkv(ji,jj) 315 311 zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) ) 316 312 zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) ) 317 ! on i-direction 318 IF ( iku > 2 ) THEN 319 ! subtract old value 320 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) 321 ! compute the new one 322 zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1) & 323 + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 324 ! add the new one to the general momentum trend 325 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) 313 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 314 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 315 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 316 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 317 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 326 318 ENDIF 327 ! on j-direction 328 IF ( ikv > 2 ) THEN 329 ! subtract old value 330 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) 331 ! compute the new one 332 zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1) & 333 + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 334 ! add the new one to the general momentum trend 335 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) 319 IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) 320 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 321 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 322 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 323 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 336 324 ENDIF 337 325 # if ! defined key_vectopt_loop … … 377 365 378 366 ! Local constant initialization 379 zcoef0 = - grav * 0.5 367 zcoef0 = - grav * 0.5_wp 380 368 ! To use density and not density anomaly 381 IF ( lk_vvl ) THEN ; znad = 1. 382 ELSE ; znad = 0. e0! Fixed volume369 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 370 ELSE ; znad = 0._wp ! Fixed volume 383 371 ENDIF 384 372 … … 392 380 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 393 381 ! s-coordinate pressure gradient correction 394 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2 *znad ) &382 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 395 383 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 396 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2 *znad ) &384 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 397 385 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 398 386 ! add to the general momentum trend … … 414 402 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 415 403 ! s-coordinate pressure gradient correction 416 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2 *znad ) &404 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 417 405 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 418 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2 *znad ) &406 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 419 407 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 420 408 ! add to the general momentum trend … … 463 451 464 452 ! Local constant initialization 465 zcoef0 = - grav * 0.5 453 zcoef0 = - grav * 0.5_wp 466 454 467 455 ! Surface value … … 495 483 & -fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) ) 496 484 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 497 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) &485 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) & 498 486 & -fse3t(ji,jj ,jk ) * rhd(ji,jj, jk) ) & 499 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &487 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) & 500 488 & -fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) ) 501 489 ! s-coordinate pressure gradient correction … … 541 529 542 530 ! Local constant initialization 543 zcoef0 = - grav * 0.5 544 zalph = 0.5 - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma545 zbeta = 0.5 + rn_gamma ! (beta =1-alpha=0.5+rn_gamma531 zcoef0 = - grav * 0.5_wp 532 zalph = 0.5_wp - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma 533 zbeta = 0.5_wp + rn_gamma ! (beta =1-alpha=0.5+rn_gamma 546 534 547 535 ! Surface value (no ponderation) … … 626 614 627 615 ! Local constant initialization 628 zcoef0 = - grav * 0.5 629 z1_10 = 1. 0 / 10.0630 z1_12 = 1. 0 / 12.0616 zcoef0 = - grav * 0.5_wp 617 z1_10 = 1._wp / 10._wp 618 z1_12 = 1._wp / 12._wp 631 619 632 620 !---------------------------------------------------------------------------------------- … … 660 648 DO jj = 2, jpjm1 661 649 DO ji = fs_2, fs_jpim1 ! vector opt. 662 cffw = 2. 0* drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1)663 664 cffu = 2. 0* drhox(ji+1,jj ,jk) * drhox(ji,jj,jk )665 cffx = 2. 0* dzx (ji+1,jj ,jk) * dzx (ji,jj,jk )650 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 651 652 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 653 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 666 654 667 cffv = 2. 0* drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk )668 cffy = 2. 0* dzy (ji ,jj+1,jk) * dzy (ji,jj,jk )655 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 656 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 669 657 670 658 IF( cffw > zep) THEN 671 drhow(ji,jj,jk) = 2. 0* drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) &672 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )659 drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 660 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 673 661 ELSE 674 drhow(ji,jj,jk) = 0. e0662 drhow(ji,jj,jk) = 0._wp 675 663 ENDIF 676 664 677 dzw(ji,jj,jk) = 2. 0* dzz(ji,jj,jk) * dzz(ji,jj,jk-1) &678 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )665 dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 666 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 679 667 680 668 IF( cffu > zep ) THEN 681 drhou(ji,jj,jk) = 2. 0* drhox(ji+1,jj,jk) * drhox(ji,jj,jk) &682 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )669 drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 670 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 683 671 ELSE 684 drhou(ji,jj,jk ) = 0. e0672 drhou(ji,jj,jk ) = 0._wp 685 673 ENDIF 686 674 687 675 IF( cffx > zep ) THEN 688 dzu(ji,jj,jk) = 2. 0*dzx(ji+1,jj,jk)*dzx(ji,jj,jk) &689 & /(dzx(ji+1,jj,jk)+dzx(ji,jj,jk))676 dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) & 677 & / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 690 678 ELSE 691 dzu(ji,jj,jk) = 0. e0679 dzu(ji,jj,jk) = 0._wp 692 680 ENDIF 693 681 694 682 IF( cffv > zep ) THEN 695 drhov(ji,jj,jk) = 2. 0* drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) &696 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )683 drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) & 684 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 697 685 ELSE 698 drhov(ji,jj,jk) = 0. e0686 drhov(ji,jj,jk) = 0._wp 699 687 ENDIF 700 688 701 689 IF( cffy > zep ) THEN 702 dzv(ji,jj,jk) = 2. 0* dzy(ji,jj+1,jk) * dzy(ji,jj,jk) &703 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )690 dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 691 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 704 692 ELSE 705 dzv(ji,jj,jk) = 0. e0693 dzv(ji,jj,jk) = 0._wp 706 694 ENDIF 707 695 … … 713 701 ! apply boundary conditions at top and bottom using 5.36-5.37 714 702 !---------------------------------------------------------------------------------- 715 drhow(:,:, 1 ) = 1.5 * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5* drhow(:,:, 2 )716 drhou(:,:, 1 ) = 1.5 * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5* drhou(:,:, 2 )717 drhov(:,:, 1 ) = 1.5 * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5* drhov(:,:, 2 )718 719 drhow(:,:,jpk) = 1.5 * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5* drhow(:,:,jpkm1)720 drhou(:,:,jpk) = 1.5 * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5* drhou(:,:,jpkm1)721 drhov(:,:,jpk) = 1.5 * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5* drhov(:,:,jpkm1)703 drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 ) 704 drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 ) 705 drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 ) 706 707 drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 708 drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 709 drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 722 710 723 711 … … 731 719 DO jj = 2, jpjm1 732 720 DO ji = fs_2, fs_jpim1 ! vector opt. 733 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &734 & * ( rhd(ji,jj,1) &735 & + 0.5 * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) &736 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) &737 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )721 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) & 722 & * ( rhd(ji,jj,1) & 723 & + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) & 724 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) & 725 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) ) 738 726 END DO 739 727 END DO … … 847 835 ! Local constant initialization 848 836 ! ------------------------------- 849 zcoef0 = - grav * 0.5 850 zforg = 0.95 e0851 zfrot = 1. e0- zforg837 zcoef0 = - grav * 0.5_wp 838 zforg = 0.95_wp 839 zfrot = 1._wp - zforg 852 840 853 841 ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance)
Note: See TracChangeset
for help on using the changeset viewer.
