- Timestamp:
- 2015-12-04T17:05:58+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5866 r6004 36 36 USE trd_oce ! trends: ocean variables 37 37 USE trddyn ! trend manager: dynamics 38 !jc USE zpshde ! partial step: hor. derivative (zps_hde routine) 38 39 ! 39 40 USE in_out_manager ! I/O manager … … 51 52 PUBLIC dyn_hpg_init ! routine called by opa module 52 53 53 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 54 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 55 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 56 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 57 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 58 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag 54 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 55 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 56 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 57 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 58 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 59 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 60 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 61 61 62 62 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) … … 131 131 !! 132 132 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 133 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf , ln_dynhpg_imp133 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 134 134 !!---------------------------------------------------------------------- 135 135 ! 136 136 REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient 137 137 READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901) 138 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp )139 138 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp ) 139 ! 140 140 REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient 141 141 READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 142 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )142 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp ) 143 143 IF(lwm) WRITE ( numond, namdyn_hpg ) 144 144 ! … … 154 154 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 155 155 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 156 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp157 156 ENDIF 158 157 ! … … 162 161 & either ln_hpg_sco or ln_hpg_prj instead') 163 162 ! 164 IF( .NOT.ln_linssh .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )&163 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 165 164 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 166 165 & ' the standard jacobian formulation hpg_sco or ' , & … … 219 218 !!---------------------------------------------------------------------- 220 219 ! 221 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )220 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 222 221 ! 223 222 IF( kt == nit000 ) THEN … … 250 249 ! hydrostatic pressure gradient 251 250 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 252 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &251 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 253 252 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 254 253 255 254 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 256 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &255 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 257 256 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 258 257 ! add to the general momentum trend … … 263 262 END DO 264 263 ! 265 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )264 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 266 265 ! 267 266 END SUBROUTINE hpg_zco … … 284 283 !!---------------------------------------------------------------------- 285 284 ! 286 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )285 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 287 286 ! 288 287 IF( kt == nit000 ) THEN … … 292 291 ENDIF 293 292 293 ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 294 !jc CALL zps_hde ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 294 295 295 296 ! Local constant initialization … … 309 310 END DO 310 311 311 312 312 ! interior value (2=<jk=<jpkm1) 313 313 DO jk = 2, jpkm1 … … 329 329 END DO 330 330 END DO 331 332 331 333 332 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) … … 353 352 END DO 354 353 ! 355 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )354 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 356 355 ! 357 356 END SUBROUTINE hpg_zps 357 358 358 359 359 SUBROUTINE hpg_sco( kt ) … … 389 389 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 390 390 ENDIF 391 392 ! Local constant initialization 391 ! 393 392 zcoef0 = - grav * 0.5_wp 394 ! To use density and not density anomaly 395 IF ( .NOT.ln_linssh ) THEN ; znad = 1._wp ! Variable volume 396 ELSE ; znad = 0._wp ! Fixed volume 393 IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 394 ELSE ; znad = 1._wp ! Variable volume: density 397 395 ENDIF 398 396 ! 399 397 ! Surface value 400 398 DO jj = 2, jpjm1 401 399 DO ji = fs_2, fs_jpim1 ! vector opt. 402 400 ! hydrostatic pressure gradient along s-surfaces 403 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) )&404 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))405 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) )&406 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))401 zhpi(ji,jj,1) = zcoef0 * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 402 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 403 zhpj(ji,jj,1) = zcoef0 * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 404 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 407 405 ! s-coordinate pressure gradient correction 408 406 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 443 441 END SUBROUTINE hpg_sco 444 442 443 445 444 SUBROUTINE hpg_isf( kt ) 446 445 !!--------------------------------------------------------------------- … … 471 470 !!---------------------------------------------------------------------- 472 471 ! 473 CALL wrk_alloc( jpi,jpj, 2, ztstop)474 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)475 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)476 ! 477 IF( kt == nit000 ) THEN472 CALL wrk_alloc( jpi,jpj, 2, ztstop ) 473 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 474 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj ) 475 ! 476 IF( kt == nit000 ) THEN 478 477 IF(lwp) WRITE(numout,*) 479 478 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 480 479 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 481 480 ENDIF 482 483 ! Local constant initialization 481 ! 484 482 zcoef0 = - grav * 0.5_wp 485 ! To use density and not density anomaly 486 ! IF ( .NOT.ln_linssh ) THEN ; znad = 1._wp ! Variable volume 487 ! ELSE ; znad = 0._wp ! Fixed volume 488 ! ENDIF 489 znad=1._wp 490 ! iniitialised to 0. zhpi zhpi 491 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 483 IF( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 484 ELSE ; znad = 1._wp ! Variable volume: density 485 ENDIF 486 zhpi(:,:,:) = 0._wp 487 zhpj(:,:,:) = 0._wp 492 488 493 489 !================================================================================== … … 496 492 497 493 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 498 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 494 ztstop(:,:,jp_tem) = -1.9_wp 495 ztstop(:,:,jp_sal) = 34.4_wp 496 497 !!gm I have the feeling that a much simplier and faster computation can be performed... 498 !!gm ====>>>> We have to discuss ! 499 500 !!gm below, faster to compute the ISF density in zrhd and remplace rhd value where tmask=0 501 !!gm furthermore, this calculation does not depends on time : do it at the first time-step only.... 499 502 500 503 ! compute density of the water displaced by the ice shelf 501 zrhd = rhd! save rhd504 zrhd(:,:,:) = rhd(:,:,:) ! save rhd 502 505 DO jk = 1, jpk 503 zdept(:,:)=gdept_1d(jk)504 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))505 END DO 506 WHERE ( tmask(:,:,:) == 1._wp)506 zdept(:,:) = gdept_1d(jk) 507 CALL eos( ztstop(:,:,:), zdept(:,:), rhd(:,:,jk) ) 508 END DO 509 WHERE( tmask(:,:,:) == 1._wp ) 507 510 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 508 511 END WHERE 509 512 510 513 ! compute rhd at the ice/oce interface (ice shelf side) 511 CALL eos( ztstop,risfdep,zrhdtop_isf)514 CALL eos( ztstop, risfdep, zrhdtop_isf ) 512 515 513 516 ! compute rhd at the ice/oce interface (ocean side) 514 DO ji =1,jpi515 DO jj =1,jpj516 ikt =mikt(ji,jj)517 ztstop(ji,jj, 1)=tsn(ji,jj,ikt,1)518 ztstop(ji,jj, 2)=tsn(ji,jj,ikt,2)517 DO ji = 1, jpi 518 DO jj = 1, jpj 519 ikt = mikt(ji,jj) 520 ztstop(ji,jj,jp_tem) = tsn(ji,jj,ikt,jp_tem) 521 ztstop(ji,jj,jp_sal) = tsn(ji,jj,ikt,jp_sal) 519 522 END DO 520 523 END DO 521 CALL eos( ztstop,risfdep,zrhdtop_oce)524 CALL eos( ztstop, risfdep, zrhdtop_oce ) 522 525 ! 523 526 ! Surface value + ice shelf gradient … … 526 529 DO jj = 1, jpj 527 530 DO ji = 1, jpi ! vector opt. 528 ikt =mikt(ji,jj)531 ikt = mikt(ji,jj) 529 532 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 530 DO jk =2,ikt-1533 DO jk = 2, ikt-1 531 534 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 532 535 & * (1._wp - tmask(ji,jj,jk)) 533 536 END DO 534 IF (ikt .GE. 2)ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) &535 &* ( risfdep(ji,jj) - gdept_1d(ikt-1) )536 END DO 537 END DO 538 riceload(:,:) = 0. 0_wp ; riceload(:,:)=ziceload(:,:)! need to be saved for diaar5537 IF( ikt >= 2 ) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 538 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 539 END DO 540 END DO 541 riceload(:,:) = 0._wp ; riceload(:,:) = ziceload(:,:) ! need to be saved for diaar5 539 542 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 540 543 DO jj = 2, jpjm1 541 544 DO ji = fs_2, fs_jpim1 ! vector opt. 542 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 545 ikt = mikt(ji,jj) 546 iktp1i = mikt(ji+1,jj) 547 iktp1j = mikt(ji,jj+1) 543 548 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 544 549 ! we assume ISF is in isostatic equilibrium 545 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj ,iktp1i) & 546 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 547 & - 0.5_wp * e3w_n(ji ,jj ,ikt ) & 548 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 549 & + ( ziceload(ji+1,jj) - ziceload(ji,jj) ) ) 550 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( 0.5_wp * e3w_n(ji ,jj+1,iktp1j) & 551 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 552 & - 0.5_wp * e3w_n(ji ,jj ,ikt ) & 553 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 554 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) 550 zhpi(ji,jj,1) = zcoef0 * ( & 551 & 0.5_wp * e3w_n(ji+1,jj,iktp1i) * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 552 & - 0.5_wp * e3w_n(ji ,jj,ikt ) * ( 2._wp * znad + rhd(ji ,jj,ikt ) + zrhdtop_oce(ji ,jj) ) & 553 & + ( ziceload(ji+1,jj) - ziceload(ji,jj) ) ) * r1_e1u(ji,jj) 554 zhpj(ji,jj,1) = zcoef0 * ( & 555 & 0.5_wp * e3w_n(ji,jj+1,iktp1j) * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 556 & - 0.5_wp * e3w_n(ji,jj ,ikt ) * ( 2._wp * znad + rhd(ji,jj ,ikt ) + zrhdtop_oce(ji,jj ) ) & 557 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) * r1_e2v(ji,jj) 555 558 ! s-coordinate pressure gradient correction (=0 if z coordinate) 556 559 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 569 572 DO ji = fs_2, fs_jpim1 ! vector opt. 570 573 iku = miku(ji,jj) 571 zpshpi(ji,jj) = 0._wp ; zpshpj(ji,jj) = 0._wp 574 zpshpi(ji,jj) = 0._wp 575 zpshpj(ji,jj) = 0._wp 572 576 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 573 577 ! u direction 574 IF ( iku .GT.1 ) THEN578 IF( iku > 1 ) THEN 575 579 ! case iku 576 zhpi(ji,jj,iku) = zcoef0 * r1_e1u(ji,jj) * ze3wu&577 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku)&578 &+ SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )580 zhpi(ji,jj,iku) = zcoef0 * r1_e1u(ji,jj) * ze3wu & 581 & * ( rhd(ji+1,jj,iku) + rhd(ji,jj,iku) & 582 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 579 583 ! corrective term ( = 0 if z coordinate ) 580 zuap 584 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj) 581 585 ! zhpi will be added in interior loop 582 ua(ji,jj,iku) 586 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap 583 587 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 584 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)= zhpi(ji,jj,iku)588 IF( mbku(ji,jj) == iku + 1 ) zpshpi(ji,jj) = zhpi(ji,jj,iku) 585 589 586 590 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 587 zhpiint = zcoef0 * r1_e1u(ji,jj)&588 & * ( e3w_n(ji+1,jj ,iku+1) * ((rhd(ji+1,jj,iku+1) + znad) &589 & 590 & - e3w_n(ji ,jj ,iku+1) * ((rhd(ji ,jj,iku+1) + znad) &591 & 592 zhpi(ji,jj,iku+1) = 591 zhpiint = zcoef0 * r1_e1u(ji,jj) & 592 & * ( e3w_n(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) & 593 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) & 594 & - e3w_n(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) & 595 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) ) 596 zhpi(ji,jj,iku+1) = zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 593 597 END IF 594 598 … … 596 600 ikv = mikv(ji,jj) 597 601 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 598 IF ( ikv .GT.1 ) THEN602 IF( ikv > 1 ) THEN 599 603 ! case ikv 600 zhpj(ji,jj,ikv) = zcoef0 * r1_e2v(ji,jj) * ze3wv&601 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv)&602 &+ SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )604 zhpj(ji,jj,ikv) = zcoef0 * r1_e2v(ji,jj) * ze3wv & 605 & * ( rhd(ji,jj+1,ikv) + rhd(ji,jj,ikv) & 606 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 603 607 ! corrective term ( = 0 if z coordinate ) 604 zvap 608 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj) 605 609 ! zhpi will be added in interior loop 606 va(ji,jj,ikv) 610 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap 607 611 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 608 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) =zhpj(ji,jj,ikv)612 IF( mbkv(ji,jj) == ikv + 1 ) zpshpj(ji,jj) = zhpj(ji,jj,ikv) 609 613 610 614 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 611 zhpjint = zcoef0 * r1_e2v(ji,jj)&612 & * ( e3w_n(ji ,jj+1,ikv+1) * ((rhd(ji,jj+1,ikv+1) + znad) &613 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv)&614 & - e3w_n(ji ,jj ,ikv+1) * ((rhd(ji,jj ,ikv+1) + znad) &615 & 615 zhpjint = zcoef0 * r1_e2v(ji,jj) & 616 & * ( e3w_n(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) & 617 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) & 618 & - e3w_n(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) & 619 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) ) 616 620 zhpj(ji,jj,ikv+1) = zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 617 END 621 ENDIF 618 622 END DO 619 623 END DO … … 969 973 ! Local constant initialization 970 974 zcoef0 = - grav 971 znad = 0.0_wp972 IF( .NOT.ln_linssh ) znad = 1._wp975 znad = 1._wp 976 IF( ln_linssh ) znad = 0._wp 973 977 974 978 ! Clean 3-D work arrays … … 1203 1207 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1204 1208 IF( .NOT.ln_linssh ) THEN 1205 1206 1209 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1210 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1207 1211 ELSE 1208 1212 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1209 1213 ENDIF 1210 1214 !!gm Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:) by definition … … 1234 1238 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1235 1239 !!---------------------------------------------------------------------- 1236 IMPLICIT NONE 1237 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 1238 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 1239 ! the interpoated function 1240 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 1241 ! 2: Linear 1240 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: fsp, xsp ! value and coordinate 1241 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1242 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1242 1243 ! 1243 1244 INTEGER :: ji, jj, jk ! dummy loop indices
Note: See TracChangeset
for help on using the changeset viewer.
