- Timestamp:
- 2017-09-27T16:29:24+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r8215 r8568 44 44 USE lib_mpp ! MPP library 45 45 USE eosbn2 ! compute density 46 USE wrk_nemo ! Memory Allocation47 46 USE timing ! Timing 48 47 USE iom … … 84 83 !!---------------------------------------------------------------------- 85 84 INTEGER, INTENT(in) :: kt ! ocean time-step index 86 REAL(wp), POINTER, DIMENSION(:,:,:) ::ztrdu, ztrdv87 !!---------------------------------------------------------------------- 88 ! 89 IF( nn_timing == 1 )CALL timing_start('dyn_hpg')85 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 86 !!---------------------------------------------------------------------- 87 ! 88 IF( ln_timing ) CALL timing_start('dyn_hpg') 90 89 ! 91 90 IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn) 92 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv)91 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 93 92 ztrdu(:,:,:) = ua(:,:,:) 94 93 ztrdv(:,:,:) = va(:,:,:) … … 108 107 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 109 108 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 110 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )109 DEALLOCATE( ztrdu , ztrdv ) 111 110 ENDIF 112 111 ! … … 114 113 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 115 114 ! 116 IF( nn_timing == 1 )CALL timing_stop('dyn_hpg')115 IF( ln_timing ) CALL timing_stop('dyn_hpg') 117 116 ! 118 117 END SUBROUTINE dyn_hpg … … 134 133 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 135 134 REAL(wp) :: znad 136 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop, zrhd! hypothesys on isf density137 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_isf! density at bottom of ISF138 REAL(wp), POINTER, DIMENSION(:,:) :: ziceload! density at bottom of ISF135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zts_top, zrhd ! hypothesys on isf density 136 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF 137 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 139 138 !! 140 139 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & … … 165 164 ! 166 165 IF( ln_hpg_djc ) & 167 & CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method & 168 & currently disabled (bugs under investigation). Please select & 169 & either ln_hpg_sco or ln_hpg_prj instead') 170 ! 171 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 172 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 173 & ' the standard jacobian formulation hpg_sco or ' , & 174 & ' the pressure jacobian formulation hpg_prj' ) 175 176 IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) & 177 & CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 178 IF( .NOT. ln_hpg_isf .AND. ln_isfcav ) & 179 & CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 166 & CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method', & 167 & ' currently disabled (bugs under investigation).' , & 168 & ' Please select either ln_hpg_sco or ln_hpg_prj instead' ) 169 ! 170 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 171 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 172 & ' the standard jacobian formulation hpg_sco or ' , & 173 & ' the pressure jacobian formulation hpg_prj' ) 174 ! 175 IF( ln_hpg_isf ) THEN 176 IF( .NOT. ln_isfcav ) CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 177 ELSE 178 IF( ln_isfcav ) CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 179 ENDIF 180 180 ! 181 181 ! ! Set nhpg from ln_hpg_... flags … … 197 197 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 198 198 ! 199 ! initialisation of ice shelf load 200 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 201 IF ( ln_isfcav ) THEN 202 CALL wrk_alloc( jpi,jpj, 2, ztstop) 203 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 204 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 199 ! 200 IF ( .NOT. ln_isfcav ) THEN !--- no ice shelf load 201 riceload(:,:) = 0._wp 202 ! 203 ELSE !--- set an ice shelf load 205 204 ! 206 205 IF(lwp) WRITE(numout,*) 207 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 208 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 209 210 ! To use density and not density anomaly 211 znad=1._wp 212 213 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 214 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 215 216 ! compute density of the water displaced by the ice shelf 217 DO jk = 1, jpk 218 CALL eos(ztstop(:,:,:),gdept_n(:,:,jk),zrhd(:,:,jk)) 219 END DO 220 221 ! compute rhd at the ice/oce interface (ice shelf side) 222 CALL eos(ztstop,risfdep,zrhdtop_isf) 223 224 ! Surface value + ice shelf gradient 225 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 226 ! divided by 2 later 227 ziceload = 0._wp 228 DO jj = 1, jpj 229 DO ji = 1, jpi 230 ikt=mikt(ji,jj) 206 IF(lwp) WRITE(numout,*) ' ice shelf case: set the ice-shelf load' 207 ALLOCATE( zts_top(jpi,jpj,jpts) , zrhd(jpi,jpj,jpk) , zrhdtop_isf(jpi,jpj) , ziceload(jpi,jpj) ) 208 ! 209 znad = 1._wp !- To use density and not density anomaly 210 ! 211 ! !- assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 212 zts_top(:,:,jp_tem) = -1.9_wp ; zts_top(:,:,jp_sal) = 34.4_wp 213 ! 214 DO jk = 1, jpk !- compute density of the water displaced by the ice shelf 215 CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) ) 216 END DO 217 ! 218 ! !- compute rhd at the ice/oce interface (ice shelf side) 219 CALL eos( zts_top , risfdep, zrhdtop_isf ) 220 ! 221 ! !- Surface value + ice shelf gradient 222 ziceload = 0._wp ! compute pressure due to ice shelf load 223 DO jj = 1, jpj ! (used to compute hpgi/j for all the level from 1 to miku/v) 224 DO ji = 1, jpi ! divided by 2 later 225 ikt = mikt(ji,jj) 231 226 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 232 DO jk =2,ikt-1227 DO jk = 2, ikt-1 233 228 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 234 229 & * (1._wp - tmask(ji,jj,jk)) 235 230 END DO 236 231 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 237 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 238 END DO 239 END DO 240 riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 241 242 CALL wrk_dealloc( jpi,jpj, 2, ztstop) 243 CALL wrk_dealloc( jpi,jpj,jpk, zrhd ) 244 CALL wrk_dealloc( jpi,jpj, zrhdtop_isf, ziceload) 245 END IF 232 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 233 END DO 234 END DO 235 riceload(:,:) = ziceload(:,:) ! need to be saved for diaar5 236 ! 237 DEALLOCATE( zts_top , zrhd , zrhdtop_isf , ziceload ) 238 ENDIF 246 239 ! 247 240 END SUBROUTINE dyn_hpg_init … … 268 261 INTEGER :: ji, jj, jk ! dummy loop indices 269 262 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars 270 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 271 !!---------------------------------------------------------------------- 272 ! 273 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 263 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 264 !!---------------------------------------------------------------------- 274 265 ! 275 266 IF( kt == nit000 ) THEN … … 315 306 END DO 316 307 ! 317 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )318 !319 308 END SUBROUTINE hpg_zco 320 309 … … 333 322 INTEGER :: iku, ikv ! temporary integers 334 323 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 335 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 336 !!---------------------------------------------------------------------- 337 ! 338 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 324 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 325 !!---------------------------------------------------------------------- 339 326 ! 340 327 IF( kt == nit000 ) THEN … … 405 392 END DO 406 393 ! 407 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )408 !409 394 END SUBROUTINE hpg_zps 410 395 … … 433 418 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 434 419 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 435 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 436 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 437 !!---------------------------------------------------------------------- 438 ! 439 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 440 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 420 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 421 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 422 !!---------------------------------------------------------------------- 441 423 ! 442 424 IF( kt == nit000 ) THEN … … 452 434 ! 453 435 IF( ln_wd ) THEN 454 DO jj = 2, jpjm1 455 DO ji = 2, jpim1 456 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 436 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 437 DO jj = 2, jpjm1 438 DO ji = 2, jpim1 439 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 440 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 458 441 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 459 442 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &443 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 461 444 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 462 445 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 463 446 464 IF(ll_tmp1) THEN465 zcpx(ji,jj) = 1.0_wp466 ELSE IF(ll_tmp2) THEN467 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here468 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj))&469 & / (sshn(ji+1,jj) - sshn(ji ,jj)))470 ELSE471 zcpx(ji,jj) = 0._wp472 ENDIF473 474 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &447 IF(ll_tmp1) THEN 448 zcpx(ji,jj) = 1.0_wp 449 ELSE IF(ll_tmp2) THEN 450 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 451 zcpx(ji,jj) = ABS( ( sshn(ji+1,jj)+ht_wd(ji+1,jj) - sshn(ji,jj)-ht_wd(ji,jj) ) & 452 & / ( sshn(ji+1,jj) - sshn(ji,jj) ) ) 453 ELSE 454 zcpx(ji,jj) = 0._wp 455 ENDIF 456 ! 457 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 475 458 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 476 459 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 477 460 & > rn_wdmin1 + rn_wdmin2 478 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &461 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 479 462 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 480 463 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 481 482 IF(ll_tmp1) THEN483 zcpy(ji,jj) = 1.0_wp484 ELSE IF(ll_tmp2) THEN485 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here486 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj))&487 & / (sshn(ji,jj+1) - sshn(ji,jj )))488 ELSE489 zcpy(ji,jj) = 0._wp490 ENDIF491 END DO492 END DO493 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )494 END 464 ! 465 IF(ll_tmp1) THEN 466 zcpy(ji,jj) = 1.0_wp 467 ELSE IF(ll_tmp2) THEN 468 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 469 zcpy(ji,jj) = ABS( ( sshn(ji,jj+1)+ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj) ) & 470 & / ( sshn(ji,jj+1) - sshn(ji,jj) ) ) 471 ELSE 472 zcpy(ji,jj) = 0._wp 473 ENDIF 474 END DO 475 END DO 476 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 477 ENDIF 495 478 496 479 ! Surface value … … 507 490 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 508 491 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 509 510 492 ! 511 493 IF( ln_wd ) THEN 512 513 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 514 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 515 zuap = zuap * zcpx(ji,jj) 516 zvap = zvap * zcpy(ji,jj) 494 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 495 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 496 zuap = zuap * zcpx(ji,jj) 497 zvap = zvap * zcpy(ji,jj) 517 498 ENDIF 518 499 ! 519 500 ! add to the general momentum trend 520 501 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 539 520 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 540 521 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 541 522 ! 542 523 IF( ln_wd ) THEN 543 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)544 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)545 zuap = zuap * zcpx(ji,jj)546 zvap = zvap * zcpy(ji,jj)547 ENDIF 548 524 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 525 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 526 zuap = zuap * zcpx(ji,jj) 527 zvap = zvap * zcpy(ji,jj) 528 ENDIF 529 ! 549 530 ! add to the general momentum trend 550 531 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 554 535 END DO 555 536 ! 556 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 557 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 537 IF( ln_wd ) DEALLOCATE( zcpx , zcpy ) 558 538 ! 559 539 END SUBROUTINE hpg_sco … … 583 563 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 584 564 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 585 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 586 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 587 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 588 !!---------------------------------------------------------------------- 589 ! 590 CALL wrk_alloc( jpi,jpj, 2, ztstop) 591 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 592 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 593 ! 594 ! Local constant initialization 595 zcoef0 = - grav * 0.5_wp 596 597 ! To use density and not density anomaly 598 znad=1._wp 599 600 ! iniitialised to 0. zhpi zhpi 601 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 565 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zhpi, zhpj 566 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts_top 567 REAL(wp), DIMENSION(jpi,jpj) :: zrhdtop_oce 568 !!---------------------------------------------------------------------- 569 ! 570 zcoef0 = - grav * 0.5_wp ! Local constant initialization 571 ! 572 znad=1._wp ! To use density and not density anomaly 573 ! 574 ! ! iniitialised to 0. zhpi zhpi 575 zhpi(:,:,:) = 0._wp ; zhpj(:,:,:) = 0._wp 602 576 603 577 ! compute rhd at the ice/oce interface (ocean side) 604 578 ! usefull to reduce residual current in the test case ISOMIP with no melting 605 DO ji =1,jpi606 DO jj =1,jpj607 ikt =mikt(ji,jj)608 zts top(ji,jj,1)=tsn(ji,jj,ikt,1)609 zts top(ji,jj,2)=tsn(ji,jj,ikt,2)579 DO ji = 1, jpi 580 DO jj = 1, jpj 581 ikt = mikt(ji,jj) 582 zts_top(ji,jj,1) = tsn(ji,jj,ikt,1) 583 zts_top(ji,jj,2) = tsn(ji,jj,ikt,2) 610 584 END DO 611 585 END DO 612 CALL eos( zts top, risfdep, zrhdtop_oce )586 CALL eos( zts_top, risfdep, zrhdtop_oce ) 613 587 614 588 !================================================================================== … … 667 641 END DO 668 642 END DO 669 !670 CALL wrk_dealloc( jpi,jpj,2 , ztstop)671 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj)672 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce )673 643 ! 674 644 END SUBROUTINE hpg_isf … … 690 660 REAL(wp) :: z1_12, cffv, cffy ! " " 691 661 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 692 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 693 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 694 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 695 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 696 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 697 !!---------------------------------------------------------------------- 698 ! 699 CALL wrk_alloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 700 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 701 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 702 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 703 ! 662 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 663 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dzx, dzy, dzz, dzu, dzv, dzw 664 REAL(wp), DIMENSION(jpi,jpj,jpk) :: drhox, drhoy, drhoz, drhou, drhov, drhow 665 REAL(wp), DIMENSION(jpi,jpj,jpk) :: rho_i, rho_j, rho_k 666 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 667 !!---------------------------------------------------------------------- 704 668 ! 705 669 IF( ln_wd ) THEN 706 DO jj = 2, jpjm1 707 DO ji = 2, jpim1 708 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 670 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 671 DO jj = 2, jpjm1 672 DO ji = 2, jpim1 673 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 709 674 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 710 675 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 711 676 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &677 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 713 678 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 714 679 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 715 680 716 IF(ll_tmp1) THEN717 zcpx(ji,jj) = 1.0_wp718 ELSE IF(ll_tmp2) THEN719 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here720 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &721 & / (sshn(ji+1,jj) - sshn(ji ,jj)) )722 ELSE723 zcpx(ji,jj) = 0._wp724 ENDIF681 IF(ll_tmp1) THEN 682 zcpx(ji,jj) = 1.0_wp 683 ELSE IF(ll_tmp2) THEN 684 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 685 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 686 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 687 ELSE 688 zcpx(ji,jj) = 0._wp 689 ENDIF 725 690 726 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &691 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 727 692 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 728 693 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 729 694 & > rn_wdmin1 + rn_wdmin2 730 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &695 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 731 696 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 732 697 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 733 698 734 IF(ll_tmp1) THEN735 zcpy(ji,jj) = 1.0_wp736 ELSE IF(ll_tmp2) THEN737 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here738 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &739 & / (sshn(ji,jj+1) - sshn(ji,jj )) )740 ELSE741 zcpy(ji,jj) = 0._wp742 ENDIF743 END DO744 END DO745 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )746 END 699 IF(ll_tmp1) THEN 700 zcpy(ji,jj) = 1.0_wp 701 ELSE IF(ll_tmp2) THEN 702 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 703 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 704 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 705 ELSE 706 zcpy(ji,jj) = 0._wp 707 ENDIF 708 END DO 709 END DO 710 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 711 ENDIF 747 712 748 713 IF( kt == nit000 ) THEN … … 903 868 END DO 904 869 END DO 905 CALL lbc_lnk( rho_k,'W',1.)906 CALL lbc_lnk( rho_i,'U',1.)907 CALL lbc_lnk( rho_j,'V',1.)870 CALL lbc_lnk( rho_k, 'W', 1. ) 871 CALL lbc_lnk( rho_i, 'U', 1. ) 872 CALL lbc_lnk( rho_j, 'V', 1. ) 908 873 909 874 … … 949 914 END DO 950 915 ! 951 CALL wrk_dealloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 952 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 953 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 954 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 916 IF( ln_wd ) DEALLOCATE( zcpx, zcpy ) 955 917 ! 956 918 END SUBROUTINE hpg_djc … … 980 942 REAL(wp) :: zrhdt1 981 943 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 982 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 983 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 984 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n 985 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 986 !!---------------------------------------------------------------------- 987 ! 988 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 989 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 990 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 991 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 944 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdept, zrhh 945 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 946 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_n, zsshv_n 947 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 948 !!---------------------------------------------------------------------- 992 949 ! 993 950 IF( kt == nit000 ) THEN … … 1003 960 1004 961 IF( ln_wd ) THEN 1005 DO jj = 2, jpjm1 1006 DO ji = 2, jpim1 1007 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 962 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 963 DO jj = 2, jpjm1 964 DO ji = 2, jpim1 965 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1008 966 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 1009 967 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 1010 968 & > rn_wdmin1 + rn_wdmin2 1011 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &969 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 1012 970 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1013 971 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1014 972 1015 IF(ll_tmp1) THEN1016 zcpx(ji,jj) = 1.0_wp1017 ELSE IF(ll_tmp2) THEN1018 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here1019 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &1020 & / (sshn(ji+1,jj) - sshn(ji ,jj)) )1021 ELSE1022 zcpx(ji,jj) = 0._wp1023 ENDIF973 IF(ll_tmp1) THEN 974 zcpx(ji,jj) = 1.0_wp 975 ELSE IF(ll_tmp2) THEN 976 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 977 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 978 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 979 ELSE 980 zcpx(ji,jj) = 0._wp 981 ENDIF 1024 982 1025 ll_tmp1 = MIN( sshn(ji,jj), sshn(ji,jj+1) ) > &983 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1026 984 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 1027 985 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 1028 986 & > rn_wdmin1 + rn_wdmin2 1029 ll_tmp2 = ( ABS( sshn(ji,jj)- sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &987 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 1030 988 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1031 989 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1032 990 1033 IF(ll_tmp1) THEN1034 zcpy(ji,jj) = 1.0_wp1035 ELSE IF(ll_tmp2) THEN1036 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here1037 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &991 IF(ll_tmp1) THEN 992 zcpy(ji,jj) = 1.0_wp 993 ELSE IF(ll_tmp2) THEN 994 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 995 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 1038 996 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1039 ELSE1040 zcpy(ji,jj) = 0._wp1041 ENDIF1042 END DO1043 END DO1044 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )1045 END 997 ELSE 998 zcpy(ji,jj) = 0._wp 999 ENDIF 1000 END DO 1001 END DO 1002 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1003 ENDIF 1046 1004 1047 1005 ! Clean 3-D work arrays … … 1298 1256 END DO 1299 1257 ! 1300 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1301 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1302 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1303 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1258 IF( ln_wd ) DEALLOCATE( zcpx, zcpy ) 1304 1259 ! 1305 1260 END SUBROUTINE hpg_prj … … 1353 1308 !!Simply geometric average 1354 1309 DO jk = 2, jpkm1-1 1355 zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1))1356 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk))1310 zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1)) 1311 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk )) 1357 1312 1358 1313 IF(zdf1 * zdf2 <= 0._wp) THEN … … 1403 1358 END DO 1404 1359 END DO 1405 1360 ! 1406 1361 ELSE 1407 1408 ENDIF 1409 1362 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1363 ENDIF 1364 ! 1410 1365 END SUBROUTINE cspline 1411 1366
Note: See TracChangeset
for help on using the changeset viewer.
