Changeset 14065
- Timestamp:
- 2020-12-03T18:48:48+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src
- Files:
-
- 10 edited
-
NST/agrif_oce_update.F90 (modified) (6 diffs)
-
OCE/DYN/dynhpg.F90 (modified) (45 diffs)
-
OCE/DYN/dynspg.F90 (modified) (7 diffs)
-
OCE/DYN/dynspg_exp.F90 (modified) (2 diffs)
-
OCE/DYN/dynspg_ts.F90 (modified) (71 diffs)
-
OCE/IOM/iom.F90 (modified) (2 diffs)
-
OCE/ISF/isf_oce.F90 (modified) (9 diffs)
-
OCE/ISF/isfload.F90 (modified) (10 diffs)
-
OCE/ISF/isfstp.F90 (modified) (9 diffs)
-
OCE/oce.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/NST/agrif_oce_update.F90
r14063 r14065 915 915 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 916 916 LOGICAL , INTENT(in ) :: before 917 !! 918 REAL(wp), DIMENSION(jpi,jpj) :: zpgu ! 2D workspace 917 919 !! 918 920 INTEGER :: ji, jj, jk … … 934 936 ! 935 937 ! Update "now" 3d velocities: 936 spgu(ji,jj) = 0._wp938 zpgu(ji,jj) = 0._wp 937 939 DO jk=1,jpkm1 938 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)940 zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 939 941 END DO 940 942 ! 941 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)943 zcorr = (tabres(ji,jj) - zpgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 942 944 DO jk=1,jpkm1 943 945 uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk) … … 954 956 ! 955 957 ! Correct "before" velocities to hold correct bt component: 956 spgu(ji,jj) = 0.e0958 zpgu(ji,jj) = 0.e0 957 959 DO jk=1,jpkm1 958 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)960 zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 959 961 END DO 960 962 ! 961 zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a)963 zcorr = uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 962 964 DO jk=1,jpkm1 963 965 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk) … … 982 984 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 983 985 LOGICAL , INTENT(in ) :: before 986 ! 987 REAL(wp), DIMENSION(jpi,jpj) :: zpgv ! 2D workspace 984 988 ! 985 989 INTEGER :: ji, jj, jk … … 1000 1004 ! 1001 1005 ! Update "now" 3d velocities: 1002 spgv(ji,jj) = 0.e01006 zpgv(ji,jj) = 0.e0 1003 1007 DO jk=1,jpkm1 1004 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)1008 zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 1005 1009 END DO 1006 1010 ! 1007 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)1011 zcorr = (tabres(ji,jj) - zpgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 1008 1012 DO jk=1,jpkm1 1009 1013 vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk) … … 1020 1024 ! 1021 1025 ! Correct "before" velocities to hold correct bt component: 1022 spgv(ji,jj) = 0.e01026 zpgv(ji,jj) = 0.e0 1023 1027 DO jk=1,jpkm1 1024 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)1028 zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 1025 1029 END DO 1026 1030 ! 1027 zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a)1031 zcorr = vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 1028 1032 DO jk=1,jpkm1 1029 1033 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk) -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynhpg.F90
r14063 r14065 44 44 USE in_out_manager ! I/O manager 45 45 USE prtctl ! Print control 46 USE lbclnk ! lateral boundary condition 46 USE lbclnk ! lateral boundary condition 47 47 USE lib_mpp ! MPP library 48 48 USE eosbn2 ! compute density … … 154 154 !! 155 155 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 156 REAL(wp) :: znad157 156 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zts_top, zrhd ! hypothesys on isf density 158 157 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF … … 206 205 ! 207 206 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 208 ! 207 ! 209 208 IF(lwp) THEN 210 209 WRITE(numout,*) … … 219 218 WRITE(numout,*) 220 219 ENDIF 221 ! 220 ! 222 221 IF ( ln_hpg_djc ) THEN 223 222 IF (ln_hpg_djc_vnh) THEN ! Von Neumann boundary condition … … 265 264 INTEGER :: ji, jj, jk ! dummy loop indices 266 265 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars 267 REAL(wp), DIMENSION(jpi,jpj ,jpk) :: zhpi, zhpj266 REAL(wp), DIMENSION(jpi,jpj) :: zhpi, zhpj 268 267 !!---------------------------------------------------------------------- 269 268 ! … … 273 272 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 274 273 ENDIF 275 276 zcoef0 = - grav * 0.5_wp ! Local constant initialization 277 278 ! Surface value 279 DO_2D( 0, 0, 0, 0 ) 274 ! 275 zcoef0 = - grav * 0.5_wp ! Local constant initialization 276 ! 277 DO_2D( 0, 0, 0, 0 ) ! Surface value 280 278 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 281 ! hydrostatic pressure gradient282 zhpi(ji,jj ,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)283 zhpj(ji,jj ,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)284 ! add to the general momentum trend285 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj ,1)286 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj ,1)279 ! ! hydrostatic pressure gradient 280 zhpi(ji,jj) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 281 zhpj(ji,jj) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 282 ! ! add to the general momentum trend 283 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj) 284 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj) 287 285 END_2D 288 289 ! 290 ! interior value (2=<jk=<jpkm1) 291 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 286 ! 287 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2=<jk=<jpkm1) 292 288 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 293 ! hydrostatic pressure gradient 294 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 295 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 296 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 297 298 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 299 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 300 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 301 ! add to the general momentum trend 302 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 303 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 289 ! ! hydrostatic pressure gradient 290 zhpi(ji,jj) = zhpi(ji,jj) + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 291 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 292 293 zhpj(ji,jj) = zhpj(ji,jj) + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 294 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 295 ! ! add to the general momentum trend 296 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj) 297 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj) 304 298 END_3D 305 299 ! … … 411 405 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 412 406 !! 413 INTEGER :: ji, jj, jk, jii, jjj 414 REAL(wp) :: zcoef0, zuap, zvap, z nad, ztmp ! temporaryscalars415 LOGICAL :: ll_tmp1, ll_tmp2 407 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 408 REAL(wp) :: zcoef0, zuap, zvap, ztmp ! local scalars 409 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 416 410 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 417 411 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter … … 427 421 ! 428 422 zcoef0 = - grav * 0.5_wp 429 IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly430 ELSE ; znad = 1._wp ! Variable volume: density431 ENDIF432 423 ! 433 424 IF( ln_wd_il ) THEN … … 450 441 zcpx(ji,jj) = 0._wp 451 442 END IF 452 443 453 444 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 454 445 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 471 462 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 472 463 END IF 473 474 ! Surface value 475 DO_2D( 0, 0, 0, 0 ) 476 ! hydrostatic pressure gradient along s-surfaces 477 zhpi(ji,jj,1) = & 478 & zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 479 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 480 & * r1_e1u(ji,jj) 481 zhpj(ji,jj,1) = & 482 & zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 483 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 484 & * r1_e2v(ji,jj) 485 ! s-coordinate pressure gradient correction 486 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 464 ! 465 DO_2D( 0, 0, 0, 0 ) ! Surface value 466 ! ! hydrostatic pressure gradient along s-surfaces 467 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) & 468 & * ( e3w(ji+1,jj ,1,Kmm) * rhd(ji+1,jj ,1) & 469 & - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) ) 470 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) & 471 & * ( e3w(ji ,jj+1,1,Kmm) * rhd(ji ,jj+1,1) & 472 & - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) ) 473 ! ! s-coordinate pressure gradient correction 474 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & 487 475 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 488 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad) &476 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & 489 477 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 490 478 ! 491 479 IF( ln_wd_il ) THEN 492 480 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 493 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 481 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 494 482 zuap = zuap * zcpx(ji,jj) 495 483 zvap = zvap * zcpy(ji,jj) 496 484 ENDIF 497 ! 498 ! add to the general momentum trend 485 ! ! add to the general momentum trend 499 486 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 500 487 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 501 488 END_2D 502 503 ! interior value (2=<jk=<jpkm1) 504 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 505 ! hydrostatic pressure gradient along s-surfaces 506 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 507 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 508 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 509 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 510 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 511 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 512 ! s-coordinate pressure gradient correction 513 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 489 ! 490 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2=<jk=<jpkm1) 491 ! ! hydrostatic pressure gradient along s-surfaces 492 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 493 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 494 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) 495 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 496 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 497 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) 498 ! ! s-coordinate pressure gradient correction 499 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) & 514 500 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 515 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad )&501 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) & 516 502 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 517 503 ! 518 504 IF( ln_wd_il ) THEN 519 505 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 520 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 506 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 521 507 zuap = zuap * zcpx(ji,jj) 522 508 zvap = zvap * zcpy(ji,jj) … … 549 535 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 550 536 !! iceload is added 551 !! 537 !! 552 538 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 553 539 !!---------------------------------------------------------------------- … … 556 542 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 557 543 !! 558 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 559 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 544 INTEGER :: ji, jj, jk ! dummy loop indices 545 INTEGER :: ikt , ikti1, iktj1 ! local integer 546 REAL(wp) :: ze3w, ze3wi1, ze3wj1 ! local scalars 547 REAL(wp) :: zcoef0, zuap, zvap ! - - 560 548 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zhpi, zhpj 561 549 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts_top … … 565 553 zcoef0 = - grav * 0.5_wp ! Local constant initialization 566 554 ! 567 znad=1._wp ! To use density and not density anomaly 568 ! 569 ! ! iniitialised to 0. zhpi zhpi 555 ! ! iniitialised to 0. zhpi zhpi 570 556 zhpi(:,:,:) = 0._wp ; zhpj(:,:,:) = 0._wp 571 557 … … 581 567 CALL eos( zts_top, risfdep, zrhdtop_oce ) 582 568 583 !================================================================================== 584 !===== Compute surface value ===================================================== 585 !================================================================================== 569 ! !===========================! 570 ! !===== surface value =====! 571 ! !===========================! 586 572 DO_2D( 0, 0, 0, 0 ) 587 ikt = mikt(ji,jj) 588 iktp1i = mikt(ji+1,jj) 589 iktp1j = mikt(ji,jj+1) 590 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 591 ! we assume ISF is in isostatic equilibrium 592 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm) & 593 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 594 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 595 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 596 & + ( risfload(ji+1,jj) - risfload(ji,jj)) ) 597 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 598 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 599 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 600 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 601 & + ( risfload(ji,jj+1) - risfload(ji,jj)) ) 602 ! s-coordinate pressure gradient correction (=0 if z coordinate) 603 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 573 ikt = mikt(ji ,jj ) ; ze3w = e3w(ji ,jj ,ikt ,Kmm) 574 ikti1 = mikt(ji+1,jj ) ; ze3wi1 = e3w(ji+1,jj ,ikti1,Kmm) 575 iktj1 = mikt(ji ,jj+1) ; ze3wj1 = e3w(ji ,jj+1,iktj1,Kmm) 576 ! ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 577 ! ! we assume ISF is in isostatic equilibrium 578 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( risfload(ji+1,jj) - risfload(ji,jj) & 579 & + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhdtop_oce(ji+1,jj) ) & 580 & - ze3w * ( rhd(ji ,jj,ikt ) + zrhdtop_oce(ji ,jj) ) ) ) 581 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( risfload(ji,jj+1) - risfload(ji,jj) & 582 & + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhdtop_oce(ji,jj+1) ) & 583 & - ze3w * ( rhd(ji,jj ,ikt ) + zrhdtop_oce(ji,jj ) ) ) ) 584 ! ! s-coordinate pressure gradient correction (=0 if z coordinate) 585 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & 604 586 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 605 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad) &587 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & 606 588 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 607 ! add to the general momentum trend589 ! ! add to the general momentum trend 608 590 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 609 591 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 610 592 END_2D 611 !================================================================================== 612 !===== Compute interior value ===================================================== 613 !================================================================================== 614 ! interior value (2=<jk=<jpkm1)593 ! 594 ! !=============================! 595 ! !===== interior values =====! 596 ! !=============================! 615 597 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 616 ! hydrostatic pressure gradient along s-surfaces 598 ze3w = e3w(ji ,jj ,jk,Kmm) 599 ze3wi1 = e3w(ji+1,jj ,jk,Kmm) 600 ze3wj1 = e3w(ji ,jj+1,jk,Kmm) 601 ! ! hydrostatic pressure gradient along s-surfaces 617 602 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 618 & * ( e3w(ji+1,jj,jk,Kmm) & 619 & * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 620 & - e3w(ji ,jj,jk,Kmm) & 621 & * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 603 & * ( ze3wi1 * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) * wmask(ji+1,jj,jk) & 604 & - ze3w * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) * wmask(ji ,jj,jk) ) 622 605 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 623 & * ( e3w(ji,jj+1,jk,Kmm) & 624 & * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 625 & - e3w(ji,jj ,jk,Kmm) & 626 & * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 627 ! s-coordinate pressure gradient correction 628 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 606 & * ( ze3wj1 * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) * wmask(ji,jj+1,jk) & 607 & - ze3w * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) * wmask(ji,jj ,jk) ) 608 ! ! s-coordinate pressure gradient correction 609 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) & 629 610 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 630 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad) &611 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) & 631 612 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 632 ! add to the general momentum trend613 ! ! add to the general momentum trend 633 614 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 634 615 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) … … 651 632 !! 652 633 INTEGER :: ji, jj, jk ! dummy loop indices 653 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 634 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 654 635 REAL(wp) :: zcoef0, zep, cffw ! temporary scalars 655 636 REAL(wp) :: z_grav_10, z1_12 … … 658 639 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 659 640 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 660 REAL(wp), DIMENSION(jpi,jpj,jpk) :: rhd_opt 661 641 662 642 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz') 663 643 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz') … … 665 645 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho') 666 646 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals 667 REAL(wp), DIMENSION(jpi,jpj) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 647 REAL(wp), DIMENSION(jpi,jpj) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 668 648 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 669 649 !!---------------------------------------------------------------------- … … 688 668 zcpx(ji,jj) = 0._wp 689 669 END IF 690 670 691 671 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 692 672 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 721 701 z1_12 = 1.0_wp / 12._wp 722 702 723 724 !! To be removed once combined with KERNEL08. KERNEL08 action removes znad from the hpg module, spg is calculated in spg moduel instead. The final code will just have rhd(:,:,:), no rhd_opt(:,:,:)725 IF( .NOT. ln_linssh ) THEN726 rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option727 ELSE728 rhd_opt(:,:,:) = rhd(:,:,:)729 END IF730 731 703 !---------------------------------------------------------------------------------------- 732 ! 1. compute and store elementary vertical differences in provisional arrays 704 ! 1. compute and store elementary vertical differences in provisional arrays 733 705 !---------------------------------------------------------------------------------------- 734 706 735 707 !!bug gm Not a true bug, but... zdzz=e3w for zdzx, zdzy verify what it is really 736 708 737 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 738 zdrhoz(ji,jj,jk) = rhd _opt (ji ,jj ,jk) - rhd_opt(ji,jj,jk-1)709 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 710 zdrhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 739 711 zdzz (ji,jj,jk) = - gde3w(ji ,jj ,jk) + gde3w(ji,jj,jk-1) 740 712 END_3D … … 745 717 zep = 1.e-15 746 718 747 !! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 719 !! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 748 720 zdrho_k(:,:,:) = 0._wp 749 721 zdz_k (:,:,:) = 0._wp 750 722 751 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 723 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 752 724 cffw = 2._wp * zdrhoz(ji ,jj ,jk) * zdrhoz(ji,jj,jk+1) 753 725 IF( cffw > zep) THEN … … 763 735 764 736 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 765 zdrho_k(:,:,1) = aco_bc_vrt * ( rhd _opt (:,:,2) - rhd_opt(:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2)737 zdrho_k(:,:,1) = aco_bc_vrt * ( rhd (:,:,2) - rhd (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 766 738 zdz_k (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k (:,:,2) 767 739 … … 769 741 IF ( mbkt(ji,jj)>1 ) THEN 770 742 iktb = mbkt(ji,jj) 771 zdrho_k(ji,jj,iktb) = aco_bc_vrt * ( rhd _opt(ji,jj,iktb) - rhd_opt(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)772 zdz_k (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k (ji,jj,iktb-1) 743 zdrho_k(ji,jj,iktb) = aco_bc_vrt * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1) 744 zdz_k (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k (ji,jj,iktb-1) 773 745 END IF 774 746 END_2D … … 778 750 !------------------------------------------------------------- 779 751 780 !! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights 781 !! rho_k stores grav * FX / rho_0 752 !! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights 753 !! rho_k stores grav * FX / rho_0 782 754 783 755 !-------------------------------------------------------------- … … 786 758 ! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1) 787 759 DO_2D( 0, 1, 0, 1) 788 z_rho_k(ji,jj,1) = grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) ) & 789 & * ( rhd _opt(ji,jj,1) &790 & + 0.5_wp * ( rhd _opt (ji,jj,2) - rhd_opt(ji,jj,1) ) &760 z_rho_k(ji,jj,1) = grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) ) & 761 & * ( rhd(ji,jj,1) & 762 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 791 763 & * ( ssh (ji,jj,Kmm) + gde3w(ji,jj,1) ) & 792 764 & / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) ) ) … … 798 770 799 771 DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 800 z_rho_k(ji,jj,jk) = zcoef0 * ( rhd _opt (ji,jj,jk) + rhd_opt(ji,jj,jk-1) ) &772 z_rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 801 773 & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) ) & 802 774 & + z_grav_10 * ( & … … 804 776 & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) & 805 777 & - ( zdz_k (ji,jj,jk) - zdz_k (ji,jj,jk-1) ) & 806 & * ( rhd _opt (ji,jj,jk) - rhd_opt(ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) &778 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) & 807 779 & ) 808 780 END_3D 809 781 810 782 !---------------------------------------------------------------------------------------- 811 ! 5. compute and store elementary horizontal differences in provisional arrays 783 ! 5. compute and store elementary horizontal differences in provisional arrays 812 784 !---------------------------------------------------------------------------------------- 813 785 814 786 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 815 zdrhox(ji,jj,jk) = rhd _opt (ji+1,jj ,jk) - rhd_opt(ji,jj,jk )787 zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 816 788 zdzx (ji,jj,jk) = - gde3w(ji+1,jj ,jk) + gde3w(ji,jj,jk ) 817 zdrhoy(ji,jj,jk) = rhd _opt (ji ,jj+1,jk) - rhd_opt(ji,jj,jk )789 zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 818 790 zdzy (ji,jj,jk) = - gde3w(ji ,jj+1,jk) + gde3w(ji,jj,jk ) 819 791 END_3D 820 792 821 CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 793 CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 822 794 823 795 !------------------------------------------------------------------------- … … 854 826 ENDIF 855 827 END_3D 856 857 !!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point 828 829 !!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point 858 830 859 831 !---------------------------------------------------------------------------------- … … 869 841 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 870 842 IF (ji < jpi) THEN 871 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 872 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd _opt (ji+1,jj,jk) - rhd_opt (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk)843 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 844 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 873 845 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) 874 846 END IF … … 877 849 IF (ji > 2) THEN 878 850 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 879 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd _opt (ji,jj,jk) - rhd_opt (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)851 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk) 880 852 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i (ji-1,jj,jk) 881 853 END IF … … 884 856 IF (jj < jpj) THEN 885 857 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 886 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd _opt (ji,jj+1,jk) - rhd_opt(ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)858 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk) 887 859 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j (ji,jj+1,jk) 888 860 END IF 889 END IF 861 END IF 890 862 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 891 863 IF (jj > 2) THEN 892 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 893 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd _opt (ji,jj,jk) - rhd_opt (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)864 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 865 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 894 866 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j (ji,jj-1,jk) 895 867 END IF … … 903 875 904 876 !-------------------------------------------------------------- 905 ! 7. Calculate integrals on side faces 877 ! 7. Calculate integrals on side faces 906 878 !------------------------------------------------------------- 907 879 908 880 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 909 881 ! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 910 z_rho_i(ji,jj,jk) = zcoef0 * ( rhd _opt (ji+1,jj,jk) + rhd_opt(ji,jj,jk) ) &911 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) 882 z_rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 883 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) 912 884 IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 913 885 z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * ( & … … 915 887 & * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) & 916 888 & - ( zdz_i (ji+1,jj,jk) - zdz_i (ji,jj,jk) ) & 917 & * ( rhd _opt (ji+1,jj,jk) - rhd_opt(ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) &889 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) & 918 890 & ) 919 891 END IF 920 921 z_rho_j(ji,jj,jk) = zcoef0 * ( rhd _opt (ji,jj+1,jk) + rhd_opt(ji,jj,jk) ) &922 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) 892 893 z_rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 894 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) 923 895 IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 924 896 z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * ( & … … 926 898 & * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) & 927 899 & - ( zdz_j (ji,jj+1,jk) - zdz_j (ji,jj,jk) ) & 928 & * ( rhd _opt (ji,jj+1,jk) - rhd_opt(ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) &900 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) & 929 901 & ) 930 902 END IF … … 932 904 933 905 !-------------------------------------------------------------- 934 ! 8. Integrate in the vertical 906 ! 8. Integrate in the vertical 935 907 !------------------------------------------------------------- 936 908 ! 937 909 ! --------------- 938 910 ! Surface value … … 943 915 IF( ln_wd_il ) THEN 944 916 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 945 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 917 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 946 918 ENDIF 947 919 ! add to the general momentum trend … … 963 935 IF( ln_wd_il ) THEN 964 936 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 965 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 937 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 966 938 ENDIF 967 939 ! add to the general momentum trend … … 1000 972 REAL(wp) :: zrhdt1 1001 973 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 974 REAL(wp), DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace 975 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_n, zsshv_n 1002 976 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdept, zrhh 1003 977 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 1004 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_n, zsshv_n1005 978 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 1006 979 !!---------------------------------------------------------------------- … … 1015 988 zcoef0 = - grav 1016 989 znad = 1._wp 1017 IF( ln_linssh ) znad = 0._wp 1018 990 IF( ln_linssh ) znad = 1._wp 991 ! 992 ! --------------- 993 ! Surface pressure gradient to be removed 994 ! --------------- 995 DO_2D( 0, 0, 0, 0 ) 996 zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 997 zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 998 END_2D 999 ! 1019 1000 IF( ln_wd_il ) THEN 1020 1001 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) … … 1034 1015 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1035 1016 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 1036 1017 1037 1018 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1038 1019 ELSE 1039 1020 zcpx(ji,jj) = 0._wp 1040 1021 END IF 1041 1022 1042 1023 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1043 1024 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 1082 1063 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 1083 1064 DO_2D( 1, 1, 1, 1 ) 1084 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad1065 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) 1085 1066 END_2D 1086 1067 … … 1120 1101 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1121 1102 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 1122 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1103 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1123 1104 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 1124 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1105 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1125 1106 !!gm not this: 1126 1107 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1127 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1108 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1128 1109 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1129 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1110 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1130 1111 END_2D 1131 1112 … … 1133 1114 1134 1115 DO_2D( 0, 0, 0, 0 ) 1135 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)1136 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad)1116 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 1117 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) ) 1137 1118 END_2D 1138 1119 … … 1216 1197 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1217 1198 ENDIF 1218 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 ) * umask(ji,jj,jk)1199 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 1219 1200 ENDIF 1220 1201 … … 1272 1253 ENDIF 1273 1254 IF( ln_wd_il ) THEN 1274 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1275 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1276 ENDIF 1277 1278 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)1255 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1256 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1257 ENDIF 1258 1259 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk) 1279 1260 ENDIF 1280 1261 ! … … 1307 1288 !!---------------------------------------------------------------------- 1308 1289 ! 1309 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1290 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1310 1291 jpi = size(fsp,1) 1311 1292 jpj = size(fsp,2) -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg.F90
r14021 r14065 82 82 INTEGER :: ji, jj, jk ! dummy loop indices 83 83 REAL(wp) :: z2dt, zg_2, zintp, zgrho0r, zld ! local scalars 84 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice 85 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 84 REAL(wp) , DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace 85 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice 86 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 86 87 !!---------------------------------------------------------------------- 87 88 ! … … 99 100 ! 100 101 DO_2D( 0, 0, 0, 0 ) 101 spgu(ji,jj) = 0._wp102 spgv(ji,jj) = 0._wp102 zpgu(ji,jj) = 0._wp 103 zpgv(ji,jj) = 0._wp 103 104 END_2D 104 105 ! … … 106 107 zg_2 = grav * 0.5 107 108 DO_2D( 0, 0, 0, 0 ) ! gradient of Patm using inverse barometer ssh 108 spgu(ji,jj) = spgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) &109 zpgu(ji,jj) = zpgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) & 109 110 & + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 110 spgv(ji,jj) = spgv(ji,jj) + zg_2 * ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) &111 zpgv(ji,jj) = zpgv(ji,jj) + zg_2 * ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) & 111 112 & + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 112 113 END_2D … … 121 122 ! 122 123 DO_2D( 0, 0, 0, 0 ) ! add tide potential forcing 123 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)124 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)124 zpgu(ji,jj) = zpgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 125 zpgv(ji,jj) = zpgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 125 126 END_2D 126 127 ! … … 128 129 zld = rn_scal_load * grav 129 130 DO_2D( 0, 0, 0, 0 ) ! add scalar approximation for load potential 130 spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)131 spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)131 zpgu(ji,jj) = zpgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 132 zpgv(ji,jj) = zpgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 132 133 END_2D 133 134 ENDIF … … 140 141 zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zgrho0r 141 142 DO_2D( 0, 0, 0, 0 ) 142 spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)143 spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)143 zpgu(ji,jj) = zpgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 144 zpgv(ji,jj) = zpgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 144 145 END_2D 145 146 DEALLOCATE( zpice ) … … 148 149 IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head 149 150 DO_2D( 0, 0, 0, 0 ) 150 spgu(ji,jj) = spgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]151 spgv(ji,jj) = spgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)151 zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3] 152 zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj) 152 153 END_2D 153 154 ENDIF 154 155 ! 155 156 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Add all terms to the general trend 156 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)157 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)157 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zpgu(ji,jj) 158 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zpgv(ji,jj) 158 159 END_3D 159 160 ! -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_exp.F90
r13497 r14065 60 60 !! 61 61 INTEGER :: ji, jj, jk ! dummy loop indices 62 REAL(wp), DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace 62 63 !!---------------------------------------------------------------------- 63 64 ! … … 67 68 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ (explicit free surface)' 68 69 ! 69 spgu(:,:) = 0._wp ; spgv(:,:) = 0._wp70 zpgu(:,:) = 0._wp ; zpgv(:,:) = 0._wp 70 71 ! 71 72 IF( .NOT.ln_linssh .AND. lwp ) WRITE(numout,*) ' non linear free surface: spg is included in dynhpg' 72 73 ENDIF 73 74 IF( ln_linssh ) THEN !* linear free surface : add the surface pressure gradient trend 75 ! 76 DO_2D( 0, 0, 0, 0 ) ! now surface pressure gradient 77 spgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 78 spgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 79 END_2D 80 ! 81 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Add it to the general trend 82 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj) 83 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj) 84 END_3D 85 ! 86 ENDIF 74 ! 75 DO_2D( 0, 0, 0, 0 ) 76 zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 77 zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 78 END_2D 79 ! 80 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 81 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zpgu(ji,jj) 82 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zpgv(ji,jj) 83 END_3D 87 84 ! 88 85 END SUBROUTINE dyn_spg_exp -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90
r14063 r14065 1 1 MODULE dynspg_ts 2 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 3 !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 4 4 5 5 !!====================================================================== … … 23 23 24 24 !!---------------------------------------------------------------------- 25 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 25 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 26 26 !! dyn_spg_ts_init: initialisation of the time-splitting scheme 27 27 !! ts_wgt : set time-splitting weights for temporal averaging (or not) … … 50 50 USE agrif_oce 51 51 #endif 52 #if defined key_asminc 52 #if defined key_asminc 53 53 USE asminc ! Assimilation increment 54 54 #endif … … 66 66 PRIVATE 67 67 68 PUBLIC dyn_spg_ts ! called by dyn_spg 68 PUBLIC dyn_spg_ts ! called by dyn_spg 69 69 PUBLIC dyn_spg_ts_init ! - - dyn_spg_init 70 70 … … 122 122 !! ** Purpose : - Compute the now trend due to the explicit time stepping 123 123 !! of the quasi-linear barotropic system, and add it to the 124 !! general momentum trend. 124 !! general momentum trend. 125 125 !! 126 126 !! ** Method : - split-explicit schem (time splitting) : 127 127 !! Barotropic variables are advanced from internal time steps 128 128 !! "n" to "n+1" if ln_bt_fw=T 129 !! or from 129 !! or from 130 130 !! "n-1" to "n+1" if ln_bt_fw=F 131 131 !! thanks to a generalized forward-backward time stepping (see ref. below). … … 136 136 !! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv 137 137 !! These are used to advect tracers and are compliant with discrete 138 !! continuity equation taken at the baroclinic time steps. This 138 !! continuity equation taken at the baroclinic time steps. This 139 139 !! ensures tracers conservation. 140 140 !! - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component. 141 141 !! 142 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 142 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 143 143 !!--------------------------------------------------------------------- 144 144 INTEGER , INTENT( in ) :: kt ! ocean time-step index … … 148 148 ! 149 149 INTEGER :: ji, jj, jk, jn ! dummy loop indices 150 LOGICAL :: ll_fw_start ! =T : forward integration 150 LOGICAL :: ll_fw_start ! =T : forward integration 151 151 LOGICAL :: ll_init ! =T : special startup of 2d equations 152 152 INTEGER :: noffset ! local integers : time offset for bdy update … … 162 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 163 163 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 164 #if defined key_qco 164 165 REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 165 ! 166 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 166 #endif 167 ! 168 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 167 169 168 170 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point … … 179 181 IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj)) 180 182 ! 181 zwdramp = r_rn_wdmin1 ! simplest ramp 183 zwdramp = r_rn_wdmin1 ! simplest ramp 182 184 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 183 ! ! inverse of baroclinic time step 184 r1_Dt_b = 1._wp / rDt 185 ! 186 ll_init = ln_bt_av ! if no time averaging, then no specific restart 185 ! ! inverse of baroclinic time step 186 r1_Dt_b = 1._wp / rDt 187 ! 188 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 189 ll_fw_start = .FALSE. 188 190 ! ! time offset in steps for bdy data update 189 191 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e 190 ELSE ; noffset = 0 192 ELSE ; noffset = 0 191 193 ENDIF 192 194 ! … … 215 217 ! and the averaging weights. We don't have an easy way of telling whether we did 216 218 ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false. 217 ! at the end of the first timestep) so just do this in all cases. 219 ! at the end of the first timestep) so just do this in all cases. 218 220 ll_fw_start = .FALSE. 219 221 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) … … 225 227 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 226 228 ! ----------------------------------------------------------------------------- 227 ! 229 ! 228 230 ! 229 231 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 230 232 ! ! --------------------------- ! 233 #if defined key_qco 231 234 DO jk = 1 , jpk 232 235 ze3u(:,:,jk) = e3u(:,:,jk,Kmm) … … 236 239 zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 237 240 zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 241 #else 242 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 243 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 244 #endif 238 245 ! 239 246 ! … … 243 250 vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk) 244 251 END DO 245 252 246 253 !!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 247 254 !!gm Is it correct to do so ? I think so... 248 249 ! != remove 2D Coriolis and pressure gradient trends=!250 ! ! -------------------------- -----------------------!255 256 ! != remove 2D Coriolis trend =! 257 ! ! -------------------------- ! 251 258 ! 252 259 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient 253 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes260 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 254 261 ! 255 262 ! !* 2D Coriolis trends 256 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 263 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 257 264 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 258 265 ! 259 266 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 260 & zu_trd, zv_trd ) ! ==>> out 261 ! 262 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) 263 ! 264 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 265 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 266 DO_2D( 0, 0, 0, 0 ) ! SPG with the application of W/D gravity filters 267 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 268 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 269 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) & 270 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 271 END_2D 272 ELSE ! now suface pressure gradient 273 DO_2D( 0, 0, 0, 0 ) 274 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) 276 END_2D 277 ENDIF 278 ! 279 ENDIF 267 & zu_trd, zv_trd ) ! ==>> out 280 268 ! 281 269 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend … … 287 275 ! ! ----------------------------------------------------------- ! 288 276 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 277 ! 289 278 ! != Add atmospheric pressure forcing =! 290 279 ! ! ---------------------------------- ! … … 319 308 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) 320 309 END_2D 321 ENDIF 310 ENDIF 322 311 ! 323 312 ! !----------------! … … 330 319 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 331 320 zztmp = r1_rho0 * r1_2 332 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:)&333 & - rnf(:,:) - rnf_b(:,:)&334 & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)&335 & + fwfisf_par(:,:) + fwfisf_par_b(:,:))321 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) & 322 & - rnf(:,:) - rnf_b(:,:) & 323 & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) & 324 & + fwfisf_par(:,:) + fwfisf_par_b(:,:) ) 336 325 ENDIF 337 326 ! != Add Stokes drift divergence =! (if exist) … … 369 358 ! 370 359 ! ----------------------------------------------------------------------- 371 ! Phase 2 : Integration of the barotropic equations 360 ! Phase 2 : Integration of the barotropic equations 372 361 ! ----------------------------------------------------------------------- 373 362 ! 374 363 ! ! ==================== ! 375 364 ! ! Initialisations ! 376 ! ! ==================== ! 377 ! Initialize barotropic variables: 365 ! ! ==================== ! 366 ! Initialize barotropic variables: 378 367 IF( ll_init )THEN 379 368 sshbb_e(:,:) = 0._wp … … 391 380 ENDIF 392 381 ! 393 IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields 394 sshn_e(:,:) = pssh (:,:,Kmm) 395 un_e (:,:) = puu_b(:,:,Kmm) 382 IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields 383 sshn_e(:,:) = pssh (:,:,Kmm) 384 un_e (:,:) = puu_b(:,:,Kmm) 396 385 vn_e (:,:) = pvv_b(:,:,Kmm) 397 386 ! 398 hu_e (:,:) = hu(:,:,Kmm) 399 hv_e (:,:) = hv(:,:,Kmm) 400 hur_e (:,:) = r1_hu(:,:,Kmm) 387 hu_e (:,:) = hu(:,:,Kmm) 388 hv_e (:,:) = hv(:,:,Kmm) 389 hur_e (:,:) = r1_hu(:,:,Kmm) 401 390 hvr_e (:,:) = r1_hv(:,:,Kmm) 402 391 ELSE ! CENTRED integration: start from BEFORE fields 403 392 sshn_e(:,:) = pssh (:,:,Kbb) 404 un_e (:,:) = puu_b(:,:,Kbb) 393 un_e (:,:) = puu_b(:,:,Kbb) 405 394 vn_e (:,:) = pvv_b(:,:,Kbb) 406 395 ! 407 hu_e (:,:) = hu(:,:,Kbb) 408 hv_e (:,:) = hv(:,:,Kbb) 409 hur_e (:,:) = r1_hu(:,:,Kbb) 396 hu_e (:,:) = hu(:,:,Kbb) 397 hv_e (:,:) = hv(:,:,Kbb) 398 hur_e (:,:) = r1_hu(:,:,Kbb) 410 399 hvr_e (:,:) = r1_hv(:,:,Kbb) 411 400 ENDIF 412 401 ! 413 402 ! Initialize sums: 414 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 403 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 415 404 pvv_b (:,:,Kaa) = 0._wp 416 405 pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level … … 419 408 ! 420 409 IF( ln_wd_dl ) THEN 421 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 422 zvwdmask(:,:) = 0._wp ! 423 zuwdav2 (:,:) = 0._wp 424 zvwdav2 (:,:) = 0._wp 425 END IF 410 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 411 zvwdmask(:,:) = 0._wp ! 412 zuwdav2 (:,:) = 0._wp 413 zvwdav2 (:,:) = 0._wp 414 END IF 426 415 427 416 ! ! ==================== ! … … 443 432 ! 444 433 ! !* Set extrapolation coefficients for predictor step: 445 IF ((jn<3).AND.ll_init) THEN ! Forward 446 za1 = 1._wp 447 za2 = 0._wp 448 za3 = 0._wp 449 ELSE ! AB3-AM4 Coefficients: bet=0.281105 434 IF ((jn<3).AND.ll_init) THEN ! Forward 435 za1 = 1._wp 436 za2 = 0._wp 437 za3 = 0._wp 438 ELSE ! AB3-AM4 Coefficients: bet=0.281105 450 439 za1 = 1.781105_wp ! za1 = 3/2 + bet 451 440 za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet) … … 467 456 !--------------------------------------------------------------------------------! 468 457 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 469 458 470 459 ! set wetting & drying mask at tracer points for this barotropic mid-step 471 460 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) … … 495 484 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 496 485 END_2D 497 #endif 486 #endif 498 487 ! 499 488 ENDIF … … 504 493 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 505 494 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 506 ! 495 ! 507 496 ! ! resulting flux at mid-step (not over the full domain) 508 497 zhU(1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column … … 515 504 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV 516 505 517 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 506 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 518 507 ! ! the direction of the flow is from dry cells 519 508 CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V 520 509 ! 521 ENDIF 510 ENDIF 522 511 ! 523 512 ! … … 543 532 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 544 533 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 545 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 534 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 546 535 IF ( ln_wd_dl_bc ) THEN 547 536 zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column … … 549 538 END IF 550 539 ! 551 ! 540 ! 552 541 ! Sea Surface Height at u-,v-points (vvl case only) 553 542 IF( .NOT.ln_linssh ) THEN … … 569 558 #endif 570 559 ENDIF 571 ! 560 ! 572 561 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 573 562 !-- m+1/2 m+1 m m-1 m-2 --! … … 626 615 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 627 616 DO_2D( 0, 0, 0, 0 ) 628 ua_e(ji,jj) = ( un_e(ji,jj) & 617 ua_e(ji,jj) = ( un_e(ji,jj) & 629 618 & + rDt_e * ( zu_spg(ji,jj) & 630 619 & + zu_trd(ji,jj) & 631 & + zu_frc(ji,jj) ) & 620 & + zu_frc(ji,jj) ) & 632 621 & ) * ssumask(ji,jj) 633 622 … … 657 646 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 658 647 ! 659 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 648 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 660 649 & + rDt_e * ( zhu_bck * zu_spg (ji,jj) & ! 661 650 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! … … 675 664 END_2D 676 665 ENDIF 677 666 678 667 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 679 668 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) … … 692 681 ! ! open boundaries 693 682 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 694 #if defined key_agrif 683 #if defined key_agrif 695 684 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif 696 685 #endif … … 711 700 ! !* Sum over whole bt loop 712 701 ! ! ---------------------- 713 za1 = wgtbtp1(jn) 702 za1 = wgtbtp1(jn) 714 703 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 715 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 716 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 704 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 705 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 717 706 ELSE ! Sum transports 718 IF ( .NOT.ln_wd_dl ) THEN 707 IF ( .NOT.ln_wd_dl ) THEN 719 708 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) 720 709 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) 721 ELSE 710 ELSE 722 711 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 723 712 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 724 END IF 713 END IF 725 714 ENDIF 726 715 ! ! Sum sea level … … 740 729 zun_save = un_adv(ji,jj) 741 730 zvn_save = vn_adv(ji,jj) 742 ! ! apply the previously computed correction 731 ! ! apply the previously computed correction 743 732 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) ) 744 733 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) ) … … 781 770 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 782 771 END_2D 783 #endif 772 #endif 784 773 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 785 774 ! … … 796 785 797 786 798 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 787 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 799 788 DO jk = 1, jpkm1 800 789 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk) … … 802 791 END DO 803 792 804 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 793 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 805 794 DO jk = 1, jpkm1 806 795 puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) & 807 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 808 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 809 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 796 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 797 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 798 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 810 799 END DO 811 END IF 812 813 800 END IF 801 802 814 803 CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current 815 804 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current … … 829 818 vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 830 819 ENDIF 831 #endif 820 #endif 832 821 ! !* write time-spliting arrays in the restart 833 822 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) … … 850 839 LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true. 851 840 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 852 INTEGER, INTENT(inout) :: jpit ! cycle length 841 INTEGER, INTENT(inout) :: jpit ! cycle length 853 842 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights 854 843 zwgt2 ! Secondary weights 855 844 856 845 INTEGER :: jic, jn, ji ! temporary integers 857 846 REAL(wp) :: za1, za2 … … 862 851 863 852 ! Set time index when averaged value is requested 864 IF (ll_fw) THEN 853 IF (ll_fw) THEN 865 854 jic = nn_e 866 855 ELSE … … 870 859 ! Set primary weights: 871 860 IF (ll_av) THEN 872 ! Define simple boxcar window for primary weights 873 ! (width = nn_e, centered around jic) 861 ! Define simple boxcar window for primary weights 862 ! (width = nn_e, centered around jic) 874 863 SELECT CASE ( nn_bt_flt ) 875 864 CASE( 0 ) ! No averaging … … 879 868 CASE( 1 ) ! Boxcar, width = nn_e 880 869 DO jn = 1, 3*nn_e 881 za1 = ABS(float(jn-jic))/float(nn_e) 870 za1 = ABS(float(jn-jic))/float(nn_e) 882 871 IF (za1 < 0.5_wp) THEN 883 872 zwgt1(jn) = 1._wp … … 888 877 CASE( 2 ) ! Boxcar, width = 2 * nn_e 889 878 DO jn = 1, 3*nn_e 890 za1 = ABS(float(jn-jic))/float(nn_e) 879 za1 = ABS(float(jn-jic))/float(nn_e) 891 880 IF (za1 < 1._wp) THEN 892 881 zwgt1(jn) = 1._wp … … 901 890 jpit = jic 902 891 ENDIF 903 892 904 893 ! Set secondary weights 905 894 DO jn = 1, jpit … … 930 919 !!---------------------------------------------------------------------- 931 920 ! 932 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 921 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 933 922 ! ! --------------- 934 923 IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN !* Read the restart file 935 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp ) 936 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp ) 937 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp ) 938 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp ) 924 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp ) 925 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp ) 926 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp ) 927 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp ) 939 928 IF( .NOT.ln_bt_av ) THEN 940 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp ) 941 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp ) 929 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp ) 930 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp ) 942 931 CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp ) 943 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp ) 944 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp ) 932 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp ) 933 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp ) 945 934 CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp ) 946 935 ENDIF … … 948 937 ! Read time integrated fluxes 949 938 IF ( .NOT.Agrif_Root() ) THEN 950 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp ) 939 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp ) 951 940 CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp ) 952 941 ELSE … … 974 963 ! 975 964 IF (.NOT.ln_bt_av) THEN 976 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 965 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 977 966 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) ) 978 967 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) ) … … 1017 1006 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1018 1007 IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax) 1019 1008 1020 1009 rDt_e = rn_Dt / REAL( nn_e , wp ) 1021 1010 zcmax = zcmax * rDt_e … … 1053 1042 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1054 1043 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e' 1055 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1044 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1056 1045 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' ) 1057 1046 END SELECT … … 1071 1060 ENDIF 1072 1061 IF( zcmax>0.9_wp ) THEN 1073 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1062 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1074 1063 ENDIF 1075 1064 ! … … 1082 1071 END SUBROUTINE dyn_spg_ts_init 1083 1072 1084 1073 1085 1074 SUBROUTINE dyn_cor_2D_init( Kmm ) 1086 1075 !!--------------------------------------------------------------------- … … 1108 1097 DO_2D( 0, 0, 0, 0 ) 1109 1098 zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) & 1110 & + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1099 & + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1111 1100 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1112 1101 END_2D … … 1145 1134 ! 1146 1135 END SELECT 1147 1136 1148 1137 END SUBROUTINE dyn_cor_2d_init 1149 1138 … … 1170 1159 ! 1171 1160 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1172 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1173 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1174 END_2D 1175 ! 1161 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1162 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1163 END_2D 1164 ! 1176 1165 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1177 1166 DO_2D( 0, 0, 0, 0 ) … … 1195 1184 END_2D 1196 1185 ! 1197 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1186 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1198 1187 DO_2D( 0, 0, 0, 0 ) 1199 1188 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & … … 1215 1204 !!---------------------------------------------------------------------- 1216 1205 !! *** ROUTINE wad_lmt *** 1217 !! 1218 !! ** Purpose : set wetting & drying mask at tracer points 1219 !! for the current barotropic sub-step 1220 !! 1221 !! ** Method : ??? 1206 !! 1207 !! ** Purpose : set wetting & drying mask at tracer points 1208 !! for the current barotropic sub-step 1209 !! 1210 !! ** Method : ??? 1222 1211 !! 1223 1212 !! ** Action : ptmsk : wetting & drying t-mask … … 1229 1218 !!---------------------------------------------------------------------- 1230 1219 ! 1231 IF( ln_wd_dl_rmp ) THEN 1220 IF( ln_wd_dl_rmp ) THEN 1232 1221 DO_2D( 1, 1, 1, 1 ) 1233 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1234 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1222 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1223 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1235 1224 ptmsk(ji,jj) = 1._wp 1236 1225 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1237 1226 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1238 ELSE 1227 ELSE 1239 1228 ptmsk(ji,jj) = 0._wp 1240 1229 ENDIF 1241 1230 END_2D 1242 ELSE 1231 ELSE 1243 1232 DO_2D( 1, 1, 1, 1 ) 1244 1233 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp … … 1254 1243 !!---------------------------------------------------------------------- 1255 1244 !! *** ROUTINE wad_lmt *** 1256 !! 1257 !! ** Purpose : set wetting & drying mask at tracer points 1258 !! for the current barotropic sub-step 1259 !! 1260 !! ** Method : ??? 1245 !! 1246 !! ** Purpose : set wetting & drying mask at tracer points 1247 !! for the current barotropic sub-step 1248 !! 1249 !! ** Method : ??? 1261 1250 !! 1262 1251 !! ** Action : ptmsk : wetting & drying t-mask … … 1270 1259 ! 1271 1260 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 1272 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1273 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1261 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1262 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1274 1263 ENDIF 1275 1264 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) … … 1279 1268 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 1280 1269 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1281 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1282 ENDIF 1283 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1270 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1271 ENDIF 1272 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1284 1273 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1285 1274 END_2D … … 1292 1281 !! *** ROUTINE wad_sp *** 1293 1282 !! 1294 !! ** Purpose : 1283 !! ** Purpose : 1295 1284 !!---------------------------------------------------------------------- 1296 1285 INTEGER :: ji ,jj ! dummy loop indices … … 1325 1314 & MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1326 1315 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1327 1316 1328 1317 IF(ll_tmp1) THEN 1329 1318 zcpy(ji,jj) = 1.0_wp … … 1337 1326 ENDIF 1338 1327 END_2D 1339 1328 1340 1329 END SUBROUTINE wad_spg 1341 1330 1342 1331 1343 1332 SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1344 1333 !!---------------------------------------------------------------------- 1345 1334 !! *** ROUTINE dyn_drg_init *** 1346 !! 1347 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1335 !! 1336 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1348 1337 !! the baroclinic part of the barotropic RHS 1349 1338 !! - compute the barotropic drag coefficients 1350 1339 !! 1351 !! ** Method : computation done over the INNER domain only 1340 !! ** Method : computation done over the INNER domain only 1352 1341 !!---------------------------------------------------------------------- 1353 1342 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices … … 1366 1355 ! 1367 1356 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities) 1368 1357 1369 1358 DO_2D( 0, 0, 0, 0 ) 1370 1359 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) … … 1381 1370 ! 1382 1371 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1383 1372 1384 1373 DO_2D( 0, 0, 0, 0 ) 1385 ikbu = mbku(ji,jj) 1386 ikbv = mbkv(ji,jj) 1374 ikbu = mbku(ji,jj) 1375 ikbv = mbkv(ji,jj) 1387 1376 zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) 1388 1377 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 1389 1378 END_2D 1390 1379 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1391 1380 1392 1381 DO_2D( 0, 0, 0, 0 ) 1393 ikbu = mbku(ji,jj) 1394 ikbv = mbkv(ji,jj) 1382 ikbu = mbku(ji,jj) 1383 ikbv = mbkv(ji,jj) 1395 1384 zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) 1396 1385 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) … … 1401 1390 zztmp = -1._wp / rDt_e 1402 1391 DO_2D( 0, 0, 0, 0 ) 1403 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1392 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1404 1393 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1405 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1394 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1406 1395 & r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1407 1396 END_2D 1408 1397 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1409 1398 1410 1399 DO_2D( 0, 0, 0, 0 ) 1411 1400 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) … … 1419 1408 ! 1420 1409 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1421 1410 1422 1411 DO_2D( 0, 0, 0, 0 ) 1423 1412 iktu = miku(ji,jj) … … 1427 1416 END_2D 1428 1417 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1429 1418 1430 1419 DO_2D( 0, 0, 0, 0 ) 1431 1420 iktu = miku(ji,jj) … … 1437 1426 ! 1438 1427 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1439 1428 1440 1429 DO_2D( 0, 0, 0, 0 ) 1441 1430 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) … … 1458 1447 ! ! set Half-step back interpolation coefficient 1459 1448 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward 1460 za0 = 1._wp 1461 za1 = 0._wp 1449 za0 = 1._wp 1450 za1 = 0._wp 1462 1451 za2 = 0._wp 1463 1452 za3 = 0._wp … … 1466 1455 za1 =-0.1666666666666_wp ! za1 = gam 1467 1456 za2 = 0.0833333333333_wp ! za2 = eps 1468 za3 = 0._wp 1469 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1470 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1457 za3 = 0._wp 1458 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1459 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1471 1460 za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps 1472 1461 za1 = 0.285_wp ! za1 = 1/2 - 2*gam - 3*eps … … 1480 1469 za2 = zgamma 1481 1470 za3 = zepsilon 1482 ENDIF 1471 ENDIF 1483 1472 ENDIF 1484 1473 END SUBROUTINE ts_bck_interp -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/IOM/iom.F90
r14063 r14065 2280 2280 ! 2281 2281 CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig0(Nis0)-1,jbegin=mjg0(Njs0)-1,ni=Ni_0,nj=Nj_0) 2282 CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj =jpj)2282 CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni=jpi, data_jbegin = -nn_hls, data_nj=jpj) 2283 2283 !don't define lon and lat for restart reading context. 2284 2284 IF ( .NOT.ldrxios ) & … … 2289 2289 ! mask land points, keep values on coast line -> specific mask for U, V and W points 2290 2290 SELECT CASE ( cdgrd ) 2291 CASE('T') ; zmask(:,:,:) = tmask(:,:,:) 2292 CASE('U') ; zmask(2:jpim1,:,:) = tmask(2:jpim1,:,:) + tmask(3:jpi,:,:) 2293 CASE('V') ; zmask(:,2:jpjm1,:) = tmask(:,2:jpjm1,:) + tmask(:,3:jpj,:) 2291 CASE('T') ; zmask( : , : ,:) = tmask(:,:,:) 2292 CASE('U') ; zmask(2:jpim1, : ,:) = tmask(2:jpim1, : ,:) + tmask(3:jpi , : ,:) 2293 CASE('V') ; zmask( : ,2:jpjm1,:) = tmask( : ,2:jpjm1,:) + tmask( : ,3:jpj,:) 2294 CASE('F') ; zmask(2:jpim1,2:jpjm1,:) = tmask(2:jpim1,2:jpjm1,:) + tmask(2:jpim1,3:jpj,:) & 2295 & + tmask(3:jpi ,2:jpjm1,:) + tmask(3:jpi ,3:jpj,:) 2294 2296 CASE('W') ; zmask(:,:,2:jpk ) = tmask(:,:,1:jpkm1) + tmask(:,:,2:jpk) ; zmask(:,:,1) = tmask(:,:,1) 2295 2297 END SELECT -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isf_oce.F90
r13558 r14065 1 1 MODULE isf_oce 2 2 !!====================================================================== 3 !! *** MODULE sbcisf***4 !! Surface module : compute iceshelf melt and heat flux3 !! *** MODULE isf_oce *** 4 !! Ice shelves : ice shelves variables defined in memory 5 5 !!====================================================================== 6 6 !! History : 3.2 ! 2011-02 (C.Harris ) Original code isf cav … … 146 146 END SUBROUTINE isf_alloc_par 147 147 148 148 149 SUBROUTINE isf_alloc_cav() 149 150 !!--------------------------------------------------------------------- … … 173 174 END SUBROUTINE isf_alloc_cav 174 175 176 175 177 SUBROUTINE isf_alloc_cpl() 176 178 !!--------------------------------------------------------------------- … … 184 186 ierr = 0 185 187 ! 186 ALLOCATE( risfcpl_ssh(jpi,jpj) , risfcpl_tsc(jpi,jpj,jpk,jpts), risfcpl_vol(jpi,jpj,jpk), STAT=ialloc )187 ierr = ierr + ialloc 188 ! 189 risfcpl_tsc(:,:,:,:) = 0. 0 ; risfcpl_vol(:,:,:) = 0.0 ; risfcpl_ssh(:,:) = 0.0190 191 IF ( ln_isfcpl_cons ) THEN192 ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) , risfcpl_cons_ssh(jpi,jpj), STAT=ialloc )188 ALLOCATE( risfcpl_ssh(jpi,jpj) , risfcpl_tsc(jpi,jpj,jpk,jpts) , risfcpl_vol(jpi,jpj,jpk) , STAT=ialloc ) 189 ierr = ierr + ialloc 190 ! 191 risfcpl_tsc(:,:,:,:) = 0._wp ; risfcpl_vol(:,:,:) = 0._wp ; risfcpl_ssh(:,:) = 0._wp 192 193 IF ( ln_isfcpl_cons ) THEN 194 ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) , risfcpl_cons_ssh(jpi,jpj) , STAT=ialloc ) 193 195 ierr = ierr + ialloc 194 196 ! 195 risfcpl_cons_tsc(:,:,:,:) = 0. 0 ; risfcpl_cons_vol(:,:,:) = 0.0 ; risfcpl_cons_ssh(:,:) = 0.0197 risfcpl_cons_tsc(:,:,:,:) = 0._wp ; risfcpl_cons_vol(:,:,:) = 0._wp ; risfcpl_cons_ssh(:,:) = 0._wp 196 198 ! 197 199 END IF … … 202 204 END SUBROUTINE isf_alloc_cpl 203 205 206 204 207 SUBROUTINE isf_dealloc_cpl() 205 208 !!--------------------------------------------------------------------- … … 213 216 ierr = 0 214 217 ! 215 DEALLOCATE( risfcpl_ssh , risfcpl_tsc, risfcpl_vol, STAT=ialloc )218 DEALLOCATE( risfcpl_ssh , risfcpl_tsc , risfcpl_vol , STAT=ialloc ) 216 219 ierr = ierr + ialloc 217 220 ! … … 221 224 END SUBROUTINE isf_dealloc_cpl 222 225 226 223 227 SUBROUTINE isf_alloc() 224 228 !!--------------------------------------------------------------------- … … 233 237 ierr = 0 ! set to zero if no array to be allocated 234 238 ! 235 ALLOCATE( fwfisf_par(jpi,jpj) , fwfisf_par_b(jpi,jpj),&236 & fwfisf_cav(jpi,jpj) , fwfisf_cav_b(jpi,jpj),&237 & fwfisf_oasis(jpi,jpj),STAT=ialloc )238 ierr = ierr + ialloc 239 ! 240 ALLOCATE( risf_par_tsc(jpi,jpj,jpts), risf_par_tsc_b(jpi,jpj,jpts), STAT=ialloc )241 ierr = ierr + ialloc 242 ! 243 ALLOCATE( risf_cav_tsc(jpi,jpj,jpts), risf_cav_tsc_b(jpi,jpj,jpts), STAT=ialloc )244 ierr = ierr + ialloc 245 ! 246 ALLOCATE( risfload(jpi,jpj), STAT=ialloc)247 ierr = ierr + ialloc 248 ! 249 ALLOCATE( mskisf_cav(jpi,jpj) , STAT=ialloc)239 ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_par_b(jpi,jpj) , & 240 & fwfisf_cav (jpi,jpj) , fwfisf_cav_b(jpi,jpj) , & 241 & fwfisf_oasis(jpi,jpj) , STAT=ialloc ) 242 ierr = ierr + ialloc 243 ! 244 ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc ) 245 ierr = ierr + ialloc 246 ! 247 ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc ) 248 ierr = ierr + ialloc 249 ! 250 ALLOCATE( risfload(jpi,jpj) , STAT=ialloc ) 251 ierr = ierr + ialloc 252 ! 253 ALLOCATE( mskisf_cav(jpi,jpj) , STAT=ialloc ) 250 254 ierr = ierr + ialloc 251 255 ! … … 254 258 ! 255 259 ! initalisation of fwf and tsc array to 0 256 risfload(:,:) = 0.0_wp 257 fwfisf_oasis(:,:) = 0.0_wp 258 fwfisf_par(:,:) = 0.0_wp ; fwfisf_par_b(:,:) = 0.0_wp 259 fwfisf_cav(:,:) = 0.0_wp ; fwfisf_cav_b(:,:) = 0.0_wp 260 risf_cav_tsc(:,:,:) = 0.0_wp ; risf_cav_tsc_b(:,:,:) = 0.0_wp 261 risf_par_tsc(:,:,:) = 0.0_wp ; risf_par_tsc_b(:,:,:) = 0.0_wp 262 ! 263 260 risfload (:,:) = 0._wp 261 fwfisf_oasis(:,:) = 0._wp 262 fwfisf_par (:,:) = 0._wp ; fwfisf_par_b (:,:) = 0._wp 263 fwfisf_cav (:,:) = 0._wp ; fwfisf_cav_b (:,:) = 0._wp 264 risf_cav_tsc(:,:,:) = 0._wp ; risf_cav_tsc_b(:,:,:) = 0._wp 265 risf_par_tsc(:,:,:) = 0._wp ; risf_par_tsc_b(:,:,:) = 0._wp 266 ! 264 267 END SUBROUTINE isf_alloc 265 268 269 !!====================================================================== 266 270 END MODULE isf_oce -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfload.F90
r13295 r14065 2 2 !!====================================================================== 3 3 !! *** MODULE isfload *** 4 !! isfload module :compute ice shelf load (needed for the hpg)4 !! Ice Shelves : compute ice shelf load (needed for the hpg) 5 5 !!====================================================================== 6 6 !! History : 4.1 ! 2019-09 (P. Mathiot) original code … … 8 8 9 9 !!---------------------------------------------------------------------- 10 !! isf load : compute ice shelf load10 !! isf_load : compute ice shelf load 11 11 !!---------------------------------------------------------------------- 12 12 … … 23 23 PRIVATE 24 24 25 PUBLIC isf_load 25 PUBLIC isf_load ! called by isfstp.F90 26 ! 26 27 !! * Substitutions 27 28 # include "do_loop_substitute.h90" 28 29 # include "domzgr_substitute.h90" 29 30 !!---------------------------------------------------------------------- 31 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 32 !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $ 33 !! Software governed by the CeCILL license (see ./LICENSE) 34 !!---------------------------------------------------------------------- 30 35 CONTAINS 31 36 … … 37 42 !! 38 43 !!-------------------------------------------------------------------- 39 !!-------------------------- OUT ------------------------------------- 40 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pisfload 41 !!-------------------------- IN ------------------------------------- 42 INTEGER, INTENT(in) :: Kmm ! ocean time level index 44 INTEGER, INTENT(in ) :: Kmm ! ocean time level index 45 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pisfload ! ice shelf load 43 46 !!---------------------------------------------------------------------- 44 47 ! … … 46 49 ! the smaller the residual flow is, the better it is. 47 50 ! 48 ! ice shelf cavity51 ! type of ice shelf cavity 49 52 SELECT CASE ( cn_isfload ) 50 53 CASE ( 'uniform' ) … … 56 59 END SUBROUTINE isf_load 57 60 58 SUBROUTINE isf_load_uniform( Kmm, pisfload ) 61 62 SUBROUTINE isf_load_uniform( Kmm, pload ) 59 63 !!-------------------------------------------------------------------- 60 64 !! *** SUBROUTINE isf_load *** … … 67 71 !! 68 72 !!-------------------------------------------------------------------- 69 !!-------------------------- OUT ------------------------------------- 70 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pisfload 71 !!-------------------------- IN ------------------------------------- 72 INTEGER, INTENT(in) :: Kmm ! ocean time level index 73 !!-------------------------------------------------------------------- 73 INTEGER, INTENT(in ) :: Kmm ! ocean time level index 74 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pload ! ice shelf load 75 ! 74 76 INTEGER :: ji, jj, jk 75 77 INTEGER :: ikt 76 REAL(wp) :: znad !77 78 REAL(wp), DIMENSION(jpi,jpj) :: zrhdtop_isf ! water density displaced by the ice shelf (at the interface) 78 79 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts_top ! water properties displaced by the ice shelf 79 80 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrhd ! water density displaced by the ice shelf 80 81 !!---------------------------------------------------------------------- 81 !82 znad = 1._wp !- To use density and not density anomaly83 82 ! 84 83 ! !- assume water displaced by the ice shelf is at T=rn_isfload_T and S=rn_isfload_S (rude) … … 87 86 DO jk = 1, jpk !- compute density of the water displaced by the ice shelf 88 87 CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) ) 88 !!st ==>> CALL eos( zts_top(:,:,:), gdept_0(:,:,jk), zrhd(:,:,jk) ) 89 89 END DO 90 90 ! … … 93 93 ! 94 94 ! !- Surface value + ice shelf gradient 95 p isfload(:,:) = 0._wp! compute pressure due to ice shelf load95 pload(:,:) = 0._wp ! compute pressure due to ice shelf load 96 96 DO_2D( 1, 1, 1, 1 ) 97 97 ikt = mikt(ji,jj) 98 98 ! 99 99 IF ( ikt > 1 ) THEN 100 ! ! top layer of the ice shelf 101 pload(ji,jj) = pload(ji,jj) & 102 & + zrhd (ji,jj,1) * e3w(ji,jj,1,Kmm) 100 103 ! 101 ! top layer of the ice shelf 102 pisfload(ji,jj) = pisfload(ji,jj) + (znad + zrhd(ji,jj,1) ) & 103 & * e3w(ji,jj,1,Kmm) 104 ! 105 ! core layers of the ice shelf 106 DO jk = 2, ikt-1 107 pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) & 108 & * e3w(ji,jj,jk,Kmm) 104 DO jk = 2, ikt-1 ! core layers of the ice shelf 105 pload(ji,jj) = pload(ji,jj) + (zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) & 106 & * e3w(ji,jj,jk,Kmm) 109 107 END DO 110 ! 111 ! deepest part of the ice shelf (between deepest T point and ice/ocean interface112 pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &113 & * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) )108 ! ! deepest part of the ice shelf (between deepest T point and ice/ocean interface 109 pload(ji,jj) = pload(ji,jj) + ( zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1) ) & 110 & * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) ) 111 !!st ==>> & * ( risfdep(ji,jj) - gdept_0(ji,jj,ikt-1) ) 114 112 ! 115 113 END IF … … 117 115 ! 118 116 END SUBROUTINE isf_load_uniform 119 117 118 !!====================================================================== 120 119 END MODULE isfload -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfstp.F90
r13237 r14065 2 2 !!====================================================================== 3 3 !! *** MODULE isfstp *** 4 !! Surface module: compute iceshelf load, melt and heat flux4 !! Ice Shelves : compute iceshelf load, melt and heat flux 5 5 !!====================================================================== 6 6 !! History : 3.2 ! 2011-02 (C.Harris ) Original code isf cav … … 42 42 !! Software governed by the CeCILL license (see ./LICENSE) 43 43 !!---------------------------------------------------------------------- 44 45 44 CONTAINS 46 45 47 SUBROUTINE isf_stp( kt, Kmm )46 SUBROUTINE isf_stp( kt, Kmm ) 48 47 !!--------------------------------------------------------------------- 49 48 !! *** ROUTINE isf_stp *** … … 58 57 !! - compute fluxes 59 58 !! - write restart variables 60 !! 61 !!---------------------------------------------------------------------- 62 INTEGER, INTENT(in) :: kt ! ocean time step 63 INTEGER, INTENT(in) :: Kmm ! ocean time level index 64 !!---------------------------------------------------------------------- 65 INTEGER :: jk ! loop index 66 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t ! e3t 59 !!---------------------------------------------------------------------- 60 INTEGER, INTENT(in) :: kt ! ocean time step 61 INTEGER, INTENT(in) :: Kmm ! ocean time level index 62 ! 63 INTEGER :: jk ! loop index 64 #if defined key_qco 65 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t ! 3D workspace 66 #endif 67 67 !!--------------------------------------------------------------------- 68 68 ! … … 83 83 ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 84 84 rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:) 85 #if defined key_qco 85 86 DO jk = 1, jpk 86 87 ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 87 88 END DO 88 CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 89 CALL isf_tbl_lvl( ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav ) 90 #else 91 CALL isf_tbl_lvl( ht(:,:), e3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav ) 92 #endif 89 93 ! 90 94 ! 1.3: compute ice shelf melt 91 CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav )95 CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav ) 92 96 ! 93 97 END IF … … 108 112 ! by simplicity, we assume the top level where param applied do not change with time (done in init part) 109 113 rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:) 114 #if defined key_qco 110 115 DO jk = 1, jpk 111 116 ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 112 117 END DO 113 CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 118 CALL isf_tbl_lvl( ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par ) 119 #else 120 CALL isf_tbl_lvl( ht(:,:), e3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par ) 121 #endif 114 122 ! 115 123 ! 2.3: compute ice shelf melt 116 CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par )124 CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par ) 117 125 ! 118 126 END IF … … 128 136 END SUBROUTINE isf_stp 129 137 130 SUBROUTINE isf_init(Kbb, Kmm, Kaa) 138 139 SUBROUTINE isf_init( Kbb, Kmm, Kaa ) 131 140 !!--------------------------------------------------------------------- 132 141 !! *** ROUTINE isfstp_init *** … … 142 151 !! - call cav/param/isfcpl init routine 143 152 !!---------------------------------------------------------------------- 144 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 153 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 154 !!---------------------------------------------------------------------- 145 155 ! 146 156 ! constrain: l_isfoasis need to be known 147 157 ! 148 ! Read namelist 149 CALL isf_nam() 150 ! 151 ! Allocate public array 152 CALL isf_alloc() 153 ! 154 ! check option compatibility 155 CALL isf_ctl() 156 ! 157 ! compute ice shelf load 158 IF ( ln_isfcav ) CALL isf_load( Kmm, risfload ) 158 CALL isf_nam() ! Read namelist 159 ! 160 CALL isf_alloc() ! Allocate public array 161 ! 162 CALL isf_ctl() ! check option compatibility 163 ! 164 IF( ln_isfcav ) CALL isf_load( Kmm, risfload ) ! compute ice shelf load 159 165 ! 160 166 ! terminate routine now if no ice shelf melt formulation specify 161 IF ( ln_isf ) THEN 162 ! 163 !--------------------------------------------------------------------------------------------------------------------- 164 ! initialisation melt in the cavity 165 IF ( ln_isfcav_mlt ) CALL isf_cav_init() 166 ! 167 !--------------------------------------------------------------------------------------------------------------------- 168 ! initialisation parametrised melt 169 IF ( ln_isfpar_mlt ) CALL isf_par_init() 170 ! 171 !--------------------------------------------------------------------------------------------------------------------- 172 ! initialisation ice sheet coupling 173 IF( ln_isfcpl ) CALL isfcpl_init(Kbb, Kmm, Kaa) 167 IF( ln_isf ) THEN 168 ! 169 IF( ln_isfcav_mlt ) CALL isf_cav_init() ! initialisation melt in the cavity 170 ! 171 IF( ln_isfpar_mlt ) CALL isf_par_init() ! initialisation parametrised melt 172 ! 173 IF( ln_isfcpl ) CALL isfcpl_init( Kbb, Kmm, Kaa ) ! initialisation ice sheet coupling 174 174 ! 175 175 END IF … … 177 177 END SUBROUTINE isf_init 178 178 179 179 180 SUBROUTINE isf_ctl() 180 181 !!--------------------------------------------------------------------- … … 283 284 END IF 284 285 END SUBROUTINE isf_ctl 285 ! 286 287 286 288 SUBROUTINE isf_nam 287 289 !!--------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/oce.F90
r14063 r14065 50 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b !: Half step time integrated fluxes 51 51 #endif 52 !53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: spgu, spgv !: horizontal surface pressure gradient54 52 55 53 !! interpolated gradient (only used in zps case) … … 91 89 ! 92 90 ierr(:) = 0 93 ALLOCATE( uu (jpi,jpj,jpk,jpt) , vv (jpi,jpj,jpk,jpt) ,&94 & ww (jpi,jpj,jpk) , hdiv(jpi,jpj,jpk) ,&95 & ts (jpi,jpj,jpk,jpts,jpt) ,&96 & rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) ,&97 & rn2b (jpi,jpj,jpk) , rn2 (jpi,jpj,jpk) ,&98 & rhd (jpi,jpj,jpk) , rhop (jpi,jpj,jpk) 91 ALLOCATE( uu (jpi,jpj,jpk,jpt) , vv (jpi,jpj,jpk,jpt) , & 92 & ww (jpi,jpj,jpk) , hdiv(jpi,jpj,jpk) , & 93 & ts (jpi,jpj,jpk,jpts,jpt) , & 94 & rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) , & 95 & rn2b (jpi,jpj,jpk) , rn2 (jpi,jpj,jpk) , & 96 & rhd (jpi,jpj,jpk) , rhop (jpi,jpj,jpk) , STAT=ierr(1) ) 99 97 ! 100 ALLOCATE( ssh(jpi,jpj,jpt) , uu_b(jpi,jpj,jpt) , vv_b(jpi,jpj,jpt) , & 101 & spgu (jpi,jpj) , spgv(jpi,jpj) , & 102 & gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts) , & 103 & gru(jpi,jpj) , grv(jpi,jpj) , & 104 & gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts) , & 105 & grui(jpi,jpj) , grvi(jpi,jpj) , & 106 & riceload(jpi,jpj) , STAT=ierr(2) ) 98 ALLOCATE( ssh (jpi,jpj,jpt) , uu_b(jpi,jpj,jpt) , vv_b(jpi,jpj,jpt) , & 99 & gtsu(jpi,jpj,jpts) , gtsv(jpi,jpj,jpts) , & 100 & gru (jpi,jpj) , grv (jpi,jpj) , & 101 & gtui(jpi,jpj,jpts) , gtvi(jpi,jpj,jpts) , & 102 & grui(jpi,jpj) , grvi(jpi,jpj) , & 103 & riceload(jpi,jpj) , STAT=ierr(2) ) 107 104 ! 108 105 ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(3) )
Note: See TracChangeset
for help on using the changeset viewer.
