Changeset 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_adv_umx.F90
- Timestamp:
- 2020-12-03T12:20:38+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_umx.F90 (modified) (51 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13292sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_adv_umx.F90
r13295 r14037 60 60 61 61 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 62 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )62 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 63 63 !!---------------------------------------------------------------------- 64 64 !! *** ROUTINE ice_dyn_adv_umx *** … … 85 85 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond concentration 86 86 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 87 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid volume 87 88 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 88 89 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 91 92 INTEGER :: icycle ! number of sub-timestep for the advection 92 93 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 93 REAL(wp) :: zdt, zvi_cen 94 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 95 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box 96 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 97 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zu_cat, zv_cat 98 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups 99 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai , z1_aip, zhvar 100 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 94 REAL(wp) :: zdt, z1_dt, zvi_cen 95 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 96 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box 97 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 98 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zu_cat, zv_cat 99 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups 100 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai , z1_aip, zhvar 101 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max 102 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max 103 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: ze_s, zes_max 101 104 ! 102 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 106 !! diagnostics 107 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 103 108 !!---------------------------------------------------------------------- 104 109 ! 105 110 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 106 111 ! 107 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 108 DO jl = 1, jpl 109 DO_2D( 0, 0, 0, 0 ) 110 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 111 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 112 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 113 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 114 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 115 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 116 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 117 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 118 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 119 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 120 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 121 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 122 END_2D 123 END DO 124 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 112 ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 113 ! thickness and salinity 114 WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 115 ELSEWHERE ; zs_i(:,:,:) = 0._wp 116 END WHERE 117 CALL icemax3D( ph_i , zhi_max ) 118 CALL icemax3D( ph_s , zhs_max ) 119 CALL icemax3D( ph_ip, zhip_max) 120 CALL icemax3D( zs_i , zsi_max ) 121 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 122 ! 123 ! enthalpies 124 DO jk = 1, nlay_i 125 WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 126 ELSEWHERE ; ze_i(:,:,jk,:) = 0._wp 127 END WHERE 128 END DO 129 DO jk = 1, nlay_s 130 WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 131 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 132 END WHERE 133 END DO 134 CALL icemax4D( ze_i , zei_max ) 135 CALL icemax4D( ze_s , zes_max ) 136 CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) 137 CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp ) 125 138 ! 126 139 ! … … 138 151 ENDIF 139 152 zdt = rDt_ice / REAL(icycle) 153 z1_dt = 1._wp / zdt 140 154 141 155 ! --- transport --- ! … … 166 180 !---------------! 167 181 DO jt = 1, icycle 182 183 ! diagnostics 184 zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 185 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 186 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 187 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 188 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 168 189 169 190 ! record at_i before advection (for open water) … … 318 339 ! 319 340 !== melt ponds ==! 320 IF ( ln_pnd_ H12) THEN341 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 321 342 ! concentration 322 343 zamsk = 1._wp … … 328 349 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 329 350 & zhvar, pv_ip, zua_ups, zva_ups ) 351 ! lid 352 IF ( ln_pnd_lids ) THEN 353 zamsk = 0._wp 354 zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:) 355 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 356 & zhvar, pv_il, zua_ups, zva_ups ) 357 ENDIF 330 358 ENDIF 359 360 ! --- Lateral boundary conditions --- ! 361 IF ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN 362 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 363 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 364 ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 365 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 366 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 367 ELSE 368 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp ) 369 ENDIF 370 CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 371 CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 331 372 ! 332 373 !== Open water area ==! … … 336 377 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 337 378 END_2D 338 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1.0_wp ) 339 ! 379 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp ) 380 ! 381 ! --- diagnostics --- ! 382 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 383 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 384 & - zdiag_adv_mass(:,:) ) * z1_dt 385 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 386 & - zdiag_adv_salt(:,:) ) * z1_dt 387 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 388 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 389 & - zdiag_adv_heat(:,:) ) * z1_dt 340 390 ! 341 391 ! --- Ensure non-negative fields and in-bound thicknesses --- ! 342 392 ! Remove negative values (conservation is ensured) 343 393 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 344 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )394 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 345 395 ! 346 396 ! --- Make sure ice thickness is not too big --- ! 347 397 ! (because ice thickness can be too large where ice concentration is very small) 348 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 398 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 399 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 349 400 ! 350 401 ! --- Ensure snow load is not too big --- ! … … 396 447 !! work on H (and not V). It is partly related to the multi-category approach 397 448 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 398 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 399 !! since sv_i and e_i are still good. 449 !! concentration is small). We also limit S and T. 400 450 !!---------------------------------------------------------------------- 401 451 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 441 491 IF( pamsk == 0._wp ) THEN 442 492 DO jl = 1, jpl 443 DO_2D( 1, 0, 1, 0 )493 DO_2D( 0, 0, 1, 0 ) 444 494 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 445 495 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) … … 450 500 ENDIF 451 501 ! 502 END_2D 503 DO_2D( 1, 0, 0, 0 ) 452 504 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 453 505 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) … … 484 536 IF( PRESENT( pua_ho ) ) THEN 485 537 DO jl = 1, jpl 486 DO_2D( 1, 0, 1, 0 ) 487 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 488 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 538 DO_2D( 0, 0, 1, 0 ) 539 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 540 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 541 END_2D 542 DO_2D( 1, 0, 0, 0 ) 543 pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 544 pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 489 545 END_2D 490 546 END DO … … 500 556 END_2D 501 557 END DO 502 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1.0_wp )503 558 ! 504 559 END SUBROUTINE adv_umx … … 539 594 ! 540 595 DO jl = 1, jpl !-- flux in x-direction 541 DO_2D( 1, 0, 1, 0 )596 DO_2D( 1, 1, 1, 0 ) 542 597 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 543 598 END_2D … … 545 600 ! 546 601 DO jl = 1, jpl !-- first guess of tracer from u-flux 547 DO_2D( 0, 0, 0, 0 )602 DO_2D( 1, 1, 0, 0 ) 548 603 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 549 604 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 552 607 END_2D 553 608 END DO 554 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )555 609 ! 556 610 DO jl = 1, jpl !-- flux in y-direction 557 DO_2D( 1, 0, 1, 0 )611 DO_2D( 1, 0, 0, 0 ) 558 612 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 559 613 END_2D … … 563 617 ! 564 618 DO jl = 1, jpl !-- flux in y-direction 565 DO_2D( 1, 0, 1, 0)619 DO_2D( 1, 0, 1, 1 ) 566 620 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 567 621 END_2D … … 569 623 ! 570 624 DO jl = 1, jpl !-- first guess of tracer from v-flux 571 DO_2D( 0, 0, 0, 0)625 DO_2D( 0, 0, 1, 1 ) 572 626 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 573 627 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 576 630 END_2D 577 631 END DO 578 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )579 632 ! 580 633 DO jl = 1, jpl !-- flux in x-direction 581 DO_2D( 1, 0, 1, 0 )634 DO_2D( 0, 0, 1, 0 ) 582 635 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 583 636 END_2D … … 628 681 ! 629 682 DO jl = 1, jpl 630 DO_2D( 1, 0, 1, 0 )683 DO_2D( 1, 1, 1, 0 ) 631 684 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 685 END_2D 686 DO_2D( 1, 0, 1, 1 ) 632 687 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 633 688 END_2D … … 646 701 ! 647 702 DO jl = 1, jpl !-- flux in x-direction 648 DO_2D( 1, 0, 1, 0 )703 DO_2D( 1, 1, 1, 0 ) 649 704 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 650 705 END_2D … … 653 708 654 709 DO jl = 1, jpl !-- first guess of tracer from u-flux 655 DO_2D( 0, 0, 0, 0 )710 DO_2D( 1, 1, 0, 0 ) 656 711 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 657 712 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 660 715 END_2D 661 716 END DO 662 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )663 717 664 718 DO jl = 1, jpl !-- flux in y-direction 665 DO_2D( 1, 0, 1, 0 )719 DO_2D( 1, 0, 0, 0 ) 666 720 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 667 721 END_2D … … 672 726 ! 673 727 DO jl = 1, jpl !-- flux in y-direction 674 DO_2D( 1, 0, 1, 0)728 DO_2D( 1, 0, 1, 1 ) 675 729 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 676 730 END_2D … … 679 733 ! 680 734 DO jl = 1, jpl !-- first guess of tracer from v-flux 681 DO_2D( 0, 0, 0, 0)735 DO_2D( 0, 0, 1, 1 ) 682 736 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 683 737 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 686 740 END_2D 687 741 END DO 688 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )689 742 ! 690 743 DO jl = 1, jpl !-- flux in x-direction 691 DO_2D( 1, 0, 1, 0 )744 DO_2D( 0, 0, 1, 0 ) 692 745 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 693 746 END_2D … … 846 899 ! 847 900 DO jl = 1, jpl 848 DO_2D( 1, 0, 1, 0 )901 DO_2D( 0, 0, 1, 0 ) 849 902 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 850 903 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 855 908 ! 856 909 DO jl = 1, jpl 857 DO_2D( 1, 0, 1, 0 )910 DO_2D( 0, 0, 1, 0 ) 858 911 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 859 912 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 865 918 ! 866 919 DO jl = 1, jpl 867 DO_2D( 1, 0, 1, 0 )920 DO_2D( 0, 0, 1, 0 ) 868 921 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 869 922 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 879 932 ! 880 933 DO jl = 1, jpl 881 DO_2D( 1, 0, 1, 0 )934 DO_2D( 0, 0, 1, 0 ) 882 935 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 883 936 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 893 946 ! 894 947 DO jl = 1, jpl 895 DO_2D( 1, 0, 1, 0 )948 DO_2D( 0, 0, 1, 0 ) 896 949 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 897 950 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 914 967 IF( ll_neg ) THEN 915 968 DO jl = 1, jpl 916 DO_2D( 1, 0, 1, 0 )969 DO_2D( 0, 0, 1, 0 ) 917 970 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 918 971 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 924 977 ! !-- High order flux in i-direction --! 925 978 DO jl = 1, jpl 926 DO_2D( 1, 0, 1, 0 )979 DO_2D( 0, 0, 1, 0 ) 927 980 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 928 981 END_2D … … 957 1010 ! !-- Laplacian in j-direction --! 958 1011 DO jl = 1, jpl 959 DO_2D( 1, 0, 0, 0 ) 1012 DO_2D( 1, 0, 0, 0 ) ! First derivative (gradient) 960 1013 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 961 1014 END_2D 962 DO_2D( 0, 0, 0, 0 ) 1015 DO_2D( 0, 0, 0, 0 ) ! Second derivative (Laplacian) 963 1016 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 964 1017 END_2D … … 968 1021 ! !-- BiLaplacian in j-direction --! 969 1022 DO jl = 1, jpl 970 DO_2D( 1, 0, 0, 0 ) 1023 DO_2D( 1, 0, 0, 0 ) ! First derivative 971 1024 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 972 1025 END_2D 973 DO_2D( 0, 0, 0, 0 ) 1026 DO_2D( 0, 0, 0, 0 ) ! Second derivative 974 1027 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 975 1028 END_2D … … 982 1035 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 983 1036 DO jl = 1, jpl 984 DO_2D( 1, 0, 1, 0 )1037 DO_2D( 1, 0, 0, 0 ) 985 1038 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 986 1039 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 990 1043 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 991 1044 DO jl = 1, jpl 992 DO_2D( 1, 0, 1, 0 )1045 DO_2D( 1, 0, 0, 0 ) 993 1046 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 994 1047 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 999 1052 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1000 1053 DO jl = 1, jpl 1001 DO_2D( 1, 0, 1, 0 )1054 DO_2D( 1, 0, 0, 0 ) 1002 1055 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1003 1056 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1012 1065 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1013 1066 DO jl = 1, jpl 1014 DO_2D( 1, 0, 1, 0 )1067 DO_2D( 1, 0, 0, 0 ) 1015 1068 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1016 1069 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1025 1078 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1026 1079 DO jl = 1, jpl 1027 DO_2D( 1, 0, 1, 0 )1080 DO_2D( 1, 0, 0, 0 ) 1028 1081 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1029 1082 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1046 1099 IF( ll_neg ) THEN 1047 1100 DO jl = 1, jpl 1048 DO_2D( 1, 0, 1, 0 )1101 DO_2D( 1, 0, 0, 0 ) 1049 1102 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1050 1103 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1056 1109 ! !-- High order flux in j-direction --! 1057 1110 DO jl = 1, jpl 1058 DO_2D( 1, 0, 1, 0 )1111 DO_2D( 1, 0, 0, 0 ) 1059 1112 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1060 1113 END_2D … … 1092 1145 ! -------------------------------------------------- 1093 1146 DO jl = 1, jpl 1094 DO_2D( 1, 0, 1, 0 )1147 DO_2D( 0, 0, 1, 0 ) 1095 1148 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1149 END_2D 1150 DO_2D( 1, 0, 0, 0 ) 1096 1151 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1097 1152 END_2D … … 1199 1254 ! --------------------------------- 1200 1255 DO jl = 1, jpl 1201 DO_2D( 1, 0, 1, 0 )1256 DO_2D( 0, 0, 1, 0 ) 1202 1257 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1203 1258 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) … … 1210 1265 END_2D 1211 1266 1212 DO_2D( 1, 0, 1, 0 )1267 DO_2D( 1, 0, 0, 0 ) 1213 1268 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1214 1269 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) … … 1409 1464 1410 1465 1411 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 1466 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 1467 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 1412 1468 !!------------------------------------------------------------------- 1413 1469 !! *** ROUTINE Hbig *** … … 1423 1479 !! ** input : Max thickness of the surrounding 9-points 1424 1480 !!------------------------------------------------------------------- 1425 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1426 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1427 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip 1481 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1482 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max, psi_max ! max ice thick from surrounding 9-pts 1483 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pes_max 1484 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pei_max 1485 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 1428 1486 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1429 ! 1430 INTEGER :: ji, jj, jl ! dummy loop indices 1431 REAL(wp) :: z1_dt, zhip, zhi, zhs, zfra 1487 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1488 ! 1489 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1490 REAL(wp) :: z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 1432 1491 !!------------------------------------------------------------------- 1433 1492 ! … … 1435 1494 ! 1436 1495 DO jl = 1, jpl 1437 1438 1496 DO_2D( 1, 1, 1, 1 ) 1439 1497 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN … … 1441 1499 ! ! -- check h_ip -- ! 1442 1500 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1443 IF( ln_pnd_H12.AND. pv_ip(ji,jj,jl) > 0._wp ) THEN1501 IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1444 1502 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1445 1503 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN … … 1468 1526 ENDIF 1469 1527 ! 1528 ! ! -- check s_i -- ! 1529 ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 1530 zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 1531 IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1532 zfra = psi_max(ji,jj,jl) / zsi 1533 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 1534 psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 1535 ENDIF 1536 ! 1470 1537 ENDIF 1471 1538 END_2D 1472 1539 END DO 1540 ! 1541 ! ! -- check e_i/v_i -- ! 1542 DO jl = 1, jpl 1543 DO_3D( 1, 1, 1, 1, 1, nlay_i ) 1544 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1545 ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 1546 zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 1547 IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1548 zfra = pei_max(ji,jj,jk,jl) / zei 1549 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1550 pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 1551 ENDIF 1552 ENDIF 1553 END_3D 1554 END DO 1555 ! ! -- check e_s/v_s -- ! 1556 DO jl = 1, jpl 1557 DO_3D( 1, 1, 1, 1, 1, nlay_s ) 1558 IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 1559 ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 1560 zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 1561 IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1562 zfra = pes_max(ji,jj,jk,jl) / zes 1563 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1564 pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 1565 ENDIF 1566 ENDIF 1567 END_3D 1568 END DO 1473 1569 ! 1474 1570 END SUBROUTINE Hbig … … 1526 1622 END SUBROUTINE Hsnow 1527 1623 1624 SUBROUTINE icemax3D( pice , pmax ) 1625 !!--------------------------------------------------------------------- 1626 !! *** ROUTINE icemax3D *** 1627 !! ** Purpose : compute the max of the 9 points around 1628 !!---------------------------------------------------------------------- 1629 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1630 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1631 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1632 INTEGER :: ji, jj, jl ! dummy loop indices 1633 !!---------------------------------------------------------------------- 1634 DO jl = 1, jpl 1635 DO jj = Njs0-1, Nje0+1 1636 DO ji = Nis0, Nie0 1637 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1638 END DO 1639 END DO 1640 DO jj = Njs0, Nje0 1641 DO ji = Nis0, Nie0 1642 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1643 END DO 1644 END DO 1645 END DO 1646 END SUBROUTINE icemax3D 1647 1648 SUBROUTINE icemax4D( pice , pmax ) 1649 !!--------------------------------------------------------------------- 1650 !! *** ROUTINE icemax4D *** 1651 !! ** Purpose : compute the max of the 9 points around 1652 !!---------------------------------------------------------------------- 1653 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1654 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1655 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1656 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1657 !!---------------------------------------------------------------------- 1658 jlay = SIZE( pice , 3 ) ! size of input arrays 1659 DO jl = 1, jpl 1660 DO jk = 1, jlay 1661 DO jj = Njs0-1, Nje0+1 1662 DO ji = Nis0, Nie0 1663 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1664 END DO 1665 END DO 1666 DO jj = Njs0, Nje0 1667 DO ji = Nis0, Nie0 1668 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1669 END DO 1670 END DO 1671 END DO 1672 END DO 1673 END SUBROUTINE icemax4D 1528 1674 1529 1675 #else
Note: See TracChangeset
for help on using the changeset viewer.
