- Timestamp:
- 2020-10-22T20:49:56+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11842_SI3-10_EAP
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_umx.F90 (modified) (48 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP
- Property svn:externals
-
old new 1 ^/utils/build/arch@HEAD arch 2 ^/utils/build/makenemo@HEAD makenemo 3 ^/utils/build/mk@HEAD mk 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 ^/vendors/FCM@HEAD ext/FCM 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 1 ^/utils/build/arch@12130 arch 2 ^/utils/build/makenemo@12191 makenemo 3 ^/utils/build/mk@11662 mk 4 ^/utils/tools_r4.0-HEAD@12672 tools 5 ^/vendors/AGRIF/dev@10586 ext/AGRIF 6 ^/vendors/FCM@10134 ext/FCM 7 ^/vendors/IOIPSL@9655 ext/IOIPSL 8 9 # SETTE mapping (inactive) 10 #^/utils/CI/sette@12135 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_adv_umx.F90
r11627 r13662 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 jj = 2, jpjm1 110 DO ji = fs_2, fs_jpim1 111 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 112 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 113 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 114 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 115 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 116 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 117 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 118 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 119 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 120 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 121 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 122 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 123 END DO 124 END DO 125 END DO 126 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 127 ! 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., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. ) 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. ) 137 CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1. ) 128 138 ! 129 139 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! … … 140 150 ENDIF 141 151 zdt = rdt_ice / REAL(icycle) 142 152 z1_dt = 1._wp / zdt 153 143 154 ! --- transport --- ! 144 155 zudy(:,:) = pu_ice(:,:) * e2u(:,:) … … 170 181 !---------------! 171 182 DO jt = 1, icycle 183 184 ! diagnostics 185 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 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 ) 172 189 173 190 ! record at_i before advection (for open water) … … 324 341 ! 325 342 !== melt ponds ==! 326 IF ( ln_pnd_ H12) THEN343 IF ( ln_pnd_LEV ) THEN 327 344 ! concentration 328 345 zamsk = 1._wp … … 334 351 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 335 352 & zhvar, pv_ip, zua_ups, zva_ups ) 353 ! lid 354 IF ( ln_pnd_lids ) THEN 355 zamsk = 0._wp 356 zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:) 357 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 358 & zhvar, pv_il, zua_ups, zva_ups ) 359 ENDIF 336 360 ENDIF 361 ! 362 ! --- Lateral boundary conditions --- ! 363 IF ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 364 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 & 365 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 366 ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 367 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 & 368 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 369 ELSE 370 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 ) 371 ENDIF 372 CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 373 CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 337 374 ! 338 375 !== Open water area ==! … … 344 381 END DO 345 382 END DO 346 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. ) 347 ! 383 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp ) 384 ! 385 ! --- diagnostics --- ! 386 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 387 & - zdiag_adv_mass(:,:) ) * z1_dt 388 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 389 & - zdiag_adv_salt(:,:) ) * z1_dt 390 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 391 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 392 & - zdiag_adv_heat(:,:) ) * z1_dt 348 393 ! 349 394 ! --- Ensure non-negative fields and in-bound thicknesses --- ! 350 395 ! Remove negative values (conservation is ensured) 351 396 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 352 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 353 ! 354 ! Make sure ice thickness is not too big 355 ! (because ice thickness can be too large where ice concentration is very small) 356 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 357 397 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 ) 398 ! 399 ! --- Make sure ice thickness is not too big --- ! 400 ! (because ice thickness can be too large where ice concentration is very small) 401 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 402 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 403 ! 404 ! --- Ensure snow load is not too big --- ! 405 CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 406 ! 358 407 END DO 359 408 ! … … 401 450 !! work on H (and not V). It is partly related to the multi-category approach 402 451 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 403 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 404 !! since sv_i and e_i are still good. 452 !! concentration is small). We also limit S and T. 405 453 !!---------------------------------------------------------------------- 406 454 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 446 494 IF( pamsk == 0._wp ) THEN 447 495 DO jl = 1, jpl 448 DO jj = 1, jpjm1496 DO jj = 2, jpjm1 449 497 DO ji = 1, fs_jpim1 450 498 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN … … 456 504 ENDIF 457 505 ! 506 END DO 507 END DO 508 DO jj = 1, jpjm1 509 DO ji = fs_2, fs_jpim1 458 510 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 459 511 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) … … 493 545 IF( PRESENT( pua_ho ) ) THEN 494 546 DO jl = 1, jpl 547 DO jj = 2, jpjm1 548 DO ji = 1, fs_jpim1 549 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 550 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 551 END DO 552 END DO 495 553 DO jj = 1, jpjm1 496 DO ji = 1, fs_jpim1497 p ua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)498 p ua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)554 DO ji = fs_2, fs_jpim1 555 pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 556 pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 499 557 END DO 500 558 END DO … … 513 571 END DO 514 572 END DO 515 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. )516 573 ! 517 574 END SUBROUTINE adv_umx … … 554 611 ! 555 612 DO jl = 1, jpl !-- flux in x-direction 556 DO jj = 1, jpj m1613 DO jj = 1, jpj 557 614 DO ji = 1, fs_jpim1 558 615 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) … … 562 619 ! 563 620 DO jl = 1, jpl !-- first guess of tracer from u-flux 564 DO jj = 2, jpjm1621 DO jj = 1, jpj 565 622 DO ji = fs_2, fs_jpim1 566 623 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & … … 571 628 END DO 572 629 END DO 573 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )574 630 ! 575 631 DO jl = 1, jpl !-- flux in y-direction 576 632 DO jj = 1, jpjm1 577 DO ji = 1, fs_jpim1633 DO ji = fs_2, fs_jpim1 578 634 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) 579 635 END DO … … 585 641 DO jl = 1, jpl !-- flux in y-direction 586 642 DO jj = 1, jpjm1 587 DO ji = 1, fs_jpim1643 DO ji = 1, jpi 588 644 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) 589 645 END DO … … 593 649 DO jl = 1, jpl !-- first guess of tracer from v-flux 594 650 DO jj = 2, jpjm1 595 DO ji = fs_2, fs_jpim1651 DO ji = 1, jpi 596 652 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 597 653 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 601 657 END DO 602 658 END DO 603 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )604 659 ! 605 660 DO jl = 1, jpl !-- flux in x-direction 606 DO jj = 1, jpjm1661 DO jj = 2, jpjm1 607 662 DO ji = 1, fs_jpim1 608 663 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) … … 657 712 ! 658 713 DO jl = 1, jpl 659 DO jj = 1, jpj m1714 DO jj = 1, jpj 660 715 DO ji = 1, fs_jpim1 661 716 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 717 END DO 718 END DO 719 DO jj = 1, jpjm1 720 DO ji = 1, jpi 662 721 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 663 722 END DO … … 677 736 ! 678 737 DO jl = 1, jpl !-- flux in x-direction 679 DO jj = 1, jpj m1738 DO jj = 1, jpj 680 739 DO ji = 1, fs_jpim1 681 740 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) … … 686 745 687 746 DO jl = 1, jpl !-- first guess of tracer from u-flux 688 DO jj = 2, jpjm1747 DO jj = 1, jpj 689 748 DO ji = fs_2, fs_jpim1 690 749 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & … … 695 754 END DO 696 755 END DO 697 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )698 756 699 757 DO jl = 1, jpl !-- flux in y-direction 700 758 DO jj = 1, jpjm1 701 DO ji = 1, fs_jpim1759 DO ji = fs_2, fs_jpim1 702 760 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 703 761 END DO … … 710 768 DO jl = 1, jpl !-- flux in y-direction 711 769 DO jj = 1, jpjm1 712 DO ji = 1, fs_jpim1770 DO ji = 1, jpi 713 771 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 714 772 END DO … … 719 777 DO jl = 1, jpl !-- first guess of tracer from v-flux 720 778 DO jj = 2, jpjm1 721 DO ji = fs_2, fs_jpim1779 DO ji = 1, jpi 722 780 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 723 781 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 727 785 END DO 728 786 END DO 729 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )730 787 ! 731 788 DO jl = 1, jpl !-- flux in x-direction 732 DO jj = 1, jpjm1789 DO jj = 2, jpjm1 733 790 DO ji = 1, fs_jpim1 734 791 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) … … 893 950 ! 894 951 DO jl = 1, jpl 895 DO jj = 1, jpjm1952 DO jj = 2, jpjm1 896 953 DO ji = 1, fs_jpim1 ! vector opt. 897 954 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 904 961 ! 905 962 DO jl = 1, jpl 906 DO jj = 1, jpjm1963 DO jj = 2, jpjm1 907 964 DO ji = 1, fs_jpim1 ! vector opt. 908 965 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) … … 916 973 ! 917 974 DO jl = 1, jpl 918 DO jj = 1, jpjm1975 DO jj = 2, jpjm1 919 976 DO ji = 1, fs_jpim1 ! vector opt. 920 977 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) … … 932 989 ! 933 990 DO jl = 1, jpl 934 DO jj = 1, jpjm1991 DO jj = 2, jpjm1 935 992 DO ji = 1, fs_jpim1 ! vector opt. 936 993 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) … … 948 1005 ! 949 1006 DO jl = 1, jpl 950 DO jj = 1, jpjm11007 DO jj = 2, jpjm1 951 1008 DO ji = 1, fs_jpim1 ! vector opt. 952 1009 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) … … 971 1028 IF( ll_neg ) THEN 972 1029 DO jl = 1, jpl 973 DO jj = 1, jpjm11030 DO jj = 2, jpjm1 974 1031 DO ji = 1, fs_jpim1 975 1032 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN … … 983 1040 ! !-- High order flux in i-direction --! 984 1041 DO jl = 1, jpl 985 DO jj = 1, jpjm11042 DO jj = 2, jpjm1 986 1043 DO ji = 1, fs_jpim1 ! vector opt. 987 1044 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) … … 1052 1109 DO jl = 1, jpl 1053 1110 DO jj = 1, jpjm1 1054 DO ji = 1, fs_jpim11111 DO ji = fs_2, fs_jpim1 1055 1112 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1056 1113 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1062 1119 DO jl = 1, jpl 1063 1120 DO jj = 1, jpjm1 1064 DO ji = 1, fs_jpim11121 DO ji = fs_2, fs_jpim1 1065 1122 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1066 1123 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 1073 1130 DO jl = 1, jpl 1074 1131 DO jj = 1, jpjm1 1075 DO ji = 1, fs_jpim11132 DO ji = fs_2, fs_jpim1 1076 1133 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1077 1134 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1088 1145 DO jl = 1, jpl 1089 1146 DO jj = 1, jpjm1 1090 DO ji = 1, fs_jpim11147 DO ji = fs_2, fs_jpim1 1091 1148 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1092 1149 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1103 1160 DO jl = 1, jpl 1104 1161 DO jj = 1, jpjm1 1105 DO ji = 1, fs_jpim11162 DO ji = fs_2, fs_jpim1 1106 1163 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1107 1164 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1126 1183 DO jl = 1, jpl 1127 1184 DO jj = 1, jpjm1 1128 DO ji = 1, fs_jpim11185 DO ji = fs_2, fs_jpim1 1129 1186 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1130 1187 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1138 1195 DO jl = 1, jpl 1139 1196 DO jj = 1, jpjm1 1140 DO ji = 1, fs_jpim1 ! vector opt.1197 DO ji = fs_2, fs_jpim1 1141 1198 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1142 1199 END DO … … 1175 1232 ! -------------------------------------------------- 1176 1233 DO jl = 1, jpl 1177 DO jj = 1, jpjm11234 DO jj = 2, jpjm1 1178 1235 DO ji = 1, fs_jpim1 ! vector opt. 1179 1236 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1237 END DO 1238 END DO 1239 DO jj = 1, jpjm1 1240 DO ji = fs_2, fs_jpim1 ! vector opt. 1180 1241 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1181 1242 END DO … … 1292 1353 ! --------------------------------- 1293 1354 DO jl = 1, jpl 1294 DO jj = 1, jpjm11355 DO jj = 2, jpjm1 1295 1356 DO ji = 1, fs_jpim1 ! vector opt. 1296 1357 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) … … 1306 1367 1307 1368 DO jj = 1, jpjm1 1308 DO ji = 1, fs_jpim1 ! vector opt.1369 DO ji = fs_2, fs_jpim1 ! vector opt. 1309 1370 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1310 1371 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) … … 1514 1575 1515 1576 1516 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 1577 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 1578 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 1517 1579 !!------------------------------------------------------------------- 1518 1580 !! *** ROUTINE Hbig *** … … 1525 1587 !! 2- check whether snow thickness is larger than the surrounding 9-points 1526 1588 !! (before advection) and reduce it by sending the excess in the ocean 1527 !! 3- check whether snow load deplets the snow-ice interface below sea level$1528 !! and reduce it by sending the excess in the ocean1529 !! 4- correct pond concentration to avoid a_ip > a_i1530 1589 !! 1531 1590 !! ** input : Max thickness of the surrounding 9-points 1532 1591 !!------------------------------------------------------------------- 1533 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1534 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1535 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 1592 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1593 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max, psi_max ! max ice thick from surrounding 9-pts 1594 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pes_max 1595 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pei_max 1596 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 1536 1597 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1537 1598 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1538 1599 ! 1539 1600 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1540 REAL(wp) :: z1_dt, zhip, zhi, zhs, zvs_excess, zfra 1541 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 1601 REAL(wp) :: z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 1542 1602 !!------------------------------------------------------------------- 1543 1603 ! … … 1545 1605 ! 1546 1606 DO jl = 1, jpl 1547 1548 1607 DO jj = 1, jpj 1549 1608 DO ji = 1, jpi … … 1552 1611 ! ! -- check h_ip -- ! 1553 1612 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1554 IF( ln_pnd_ H12.AND. pv_ip(ji,jj,jl) > 0._wp ) THEN1613 IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1555 1614 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1556 1615 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN … … 1578 1637 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1579 1638 ENDIF 1639 ! 1640 ! ! -- check s_i -- ! 1641 ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 1642 zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 1643 IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1644 zfra = psi_max(ji,jj,jl) / zsi 1645 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 1646 psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 1647 ENDIF 1580 1648 ! 1581 ! ! -- check snow load -- ! 1582 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 1583 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 1584 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 1649 ENDIF 1650 END DO 1651 END DO 1652 END DO 1653 ! 1654 ! ! -- check e_i/v_i -- ! 1655 DO jl = 1, jpl 1656 DO jk = 1, nlay_i 1657 DO jj = 1, jpj 1658 DO ji = 1, jpi 1659 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1660 ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 1661 zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 1662 IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1663 zfra = pei_max(ji,jj,jk,jl) / zei 1664 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1665 pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 1666 ENDIF 1667 ENDIF 1668 END DO 1669 END DO 1670 END DO 1671 END DO 1672 ! ! -- check e_s/v_s -- ! 1673 DO jl = 1, jpl 1674 DO jk = 1, nlay_s 1675 DO jj = 1, jpj 1676 DO ji = 1, jpi 1677 IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 1678 ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 1679 zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 1680 IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1681 zfra = pes_max(ji,jj,jk,jl) / zes 1682 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1683 pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 1684 ENDIF 1685 ENDIF 1686 END DO 1687 END DO 1688 END DO 1689 END DO 1690 ! 1691 END SUBROUTINE Hbig 1692 1693 1694 SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 1695 !!------------------------------------------------------------------- 1696 !! *** ROUTINE Hsnow *** 1697 !! 1698 !! ** Purpose : 1- Check snow load after advection 1699 !! 2- Correct pond concentration to avoid a_ip > a_i 1700 !! 1701 !! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface 1702 !! then put the snow excess in the ocean 1703 !! 1704 !! ** Notes : This correction is crucial because of the call to routine icecor afterwards 1705 !! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 1706 !! make the snow very thick (if concentration decreases drastically) 1707 !! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 1708 !!------------------------------------------------------------------- 1709 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1710 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip 1711 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1712 ! 1713 INTEGER :: ji, jj, jl ! dummy loop indices 1714 REAL(wp) :: z1_dt, zvs_excess, zfra 1715 !!------------------------------------------------------------------- 1716 ! 1717 z1_dt = 1._wp / pdt 1718 ! 1719 ! -- check snow load -- ! 1720 DO jl = 1, jpl 1721 DO jj = 1, jpj 1722 DO ji = 1, jpi 1723 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1724 ! 1585 1725 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 1586 IF( zvs_excess > 0._wp ) THEN 1726 ! 1727 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 1728 ! put snow excess in the ocean 1587 1729 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1588 1730 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1589 1731 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1590 ! 1732 ! correct snow volume and heat content 1591 1733 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1592 1734 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1593 1735 ENDIF 1594 1736 ! 1595 1737 ENDIF 1596 1738 END DO 1597 1739 END DO 1598 END DO 1599 ! !-- correct pond concentration to avoid a_ip > a_i 1740 END DO 1741 ! 1742 !-- correct pond concentration to avoid a_ip > a_i -- ! 1600 1743 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1601 1744 ! 1602 ! 1603 END SUBROUTINE Hbig 1604 1745 END SUBROUTINE Hsnow 1746 1747 SUBROUTINE icemax3D( pice , pmax ) 1748 !!--------------------------------------------------------------------- 1749 !! *** ROUTINE icemax3D *** 1750 !! ** Purpose : compute the max of the 9 points around 1751 !!---------------------------------------------------------------------- 1752 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1753 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1754 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1755 INTEGER :: ji, jj, jl ! dummy loop indices 1756 !!---------------------------------------------------------------------- 1757 DO jl = 1, jpl 1758 DO jj = 1, jpj 1759 DO ji = 2, jpim1 1760 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1761 END DO 1762 END DO 1763 DO jj = 2, jpjm1 1764 DO ji = 2, jpim1 1765 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1766 END DO 1767 END DO 1768 END DO 1769 END SUBROUTINE icemax3D 1770 1771 SUBROUTINE icemax4D( pice , pmax ) 1772 !!--------------------------------------------------------------------- 1773 !! *** ROUTINE icemax4D *** 1774 !! ** Purpose : compute the max of the 9 points around 1775 !!---------------------------------------------------------------------- 1776 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1777 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1778 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1779 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1780 !!---------------------------------------------------------------------- 1781 jlay = SIZE( pice , 3 ) ! size of input arrays 1782 DO jl = 1, jpl 1783 DO jk = 1, jlay 1784 DO jj = 1, jpj 1785 DO ji = 2, jpim1 1786 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1787 END DO 1788 END DO 1789 DO jj = 2, jpjm1 1790 DO ji = 2, jpim1 1791 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1792 END DO 1793 END DO 1794 END DO 1795 END DO 1796 END SUBROUTINE icemax4D 1797 1605 1798 #else 1606 1799 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
