- Timestamp:
- 2020-12-02T14:55:21+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_umx.F90 (modified) (51 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette@13 292sette10 ^/utils/CI/sette@13559 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/ICE/icedyn_adv_umx.F90
r13295 r13998 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 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 186 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 187 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 168 188 169 189 ! record at_i before advection (for open water) … … 318 338 ! 319 339 !== melt ponds ==! 320 IF ( ln_pnd_ H12) THEN340 IF ( ln_pnd_LEV ) THEN 321 341 ! concentration 322 342 zamsk = 1._wp … … 328 348 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 329 349 & zhvar, pv_ip, zua_ups, zva_ups ) 350 ! lid 351 IF ( ln_pnd_lids ) THEN 352 zamsk = 0._wp 353 zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:) 354 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 355 & zhvar, pv_il, zua_ups, zva_ups ) 356 ENDIF 330 357 ENDIF 358 359 ! --- Lateral boundary conditions --- ! 360 IF ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 361 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 & 362 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 363 ELSEIF( ln_pnd_LEV .AND. .NOT.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 ) 366 ELSE 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 ENDIF 369 CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 370 CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 331 371 ! 332 372 !== Open water area ==! … … 336 376 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 337 377 END_2D 338 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1.0_wp ) 339 ! 378 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp ) 379 ! 380 ! --- diagnostics --- ! 381 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 382 & - zdiag_adv_mass(:,:) ) * z1_dt 383 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 384 & - zdiag_adv_salt(:,:) ) * z1_dt 385 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 386 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 387 & - zdiag_adv_heat(:,:) ) * z1_dt 340 388 ! 341 389 ! --- Ensure non-negative fields and in-bound thicknesses --- ! 342 390 ! Remove negative values (conservation is ensured) 343 391 ! (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 )392 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 393 ! 346 394 ! --- Make sure ice thickness is not too big --- ! 347 395 ! (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 ) 396 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 397 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 349 398 ! 350 399 ! --- Ensure snow load is not too big --- ! … … 396 445 !! work on H (and not V). It is partly related to the multi-category approach 397 446 !! 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. 447 !! concentration is small). We also limit S and T. 400 448 !!---------------------------------------------------------------------- 401 449 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 441 489 IF( pamsk == 0._wp ) THEN 442 490 DO jl = 1, jpl 443 DO_2D( 1, 0, 1, 0 )491 DO_2D( 0, 0, 1, 0 ) 444 492 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 445 493 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) … … 450 498 ENDIF 451 499 ! 500 END_2D 501 DO_2D( 1, 0, 0, 0 ) 452 502 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 453 503 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) … … 484 534 IF( PRESENT( pua_ho ) ) THEN 485 535 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) 536 DO_2D( 0, 0, 1, 0 ) 537 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 538 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 539 END_2D 540 DO_2D( 1, 0, 0, 0 ) 541 pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 542 pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 489 543 END_2D 490 544 END DO … … 500 554 END_2D 501 555 END DO 502 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1.0_wp )503 556 ! 504 557 END SUBROUTINE adv_umx … … 539 592 ! 540 593 DO jl = 1, jpl !-- flux in x-direction 541 DO_2D( 1, 0, 1, 0 )594 DO_2D( 1, 1, 1, 0 ) 542 595 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 596 END_2D … … 545 598 ! 546 599 DO jl = 1, jpl !-- first guess of tracer from u-flux 547 DO_2D( 0, 0, 0, 0 )600 DO_2D( 1, 1, 0, 0 ) 548 601 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 549 602 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 552 605 END_2D 553 606 END DO 554 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )555 607 ! 556 608 DO jl = 1, jpl !-- flux in y-direction 557 DO_2D( 1, 0, 1, 0 )609 DO_2D( 1, 0, 0, 0 ) 558 610 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 611 END_2D … … 563 615 ! 564 616 DO jl = 1, jpl !-- flux in y-direction 565 DO_2D( 1, 0, 1, 0)617 DO_2D( 1, 0, 1, 1 ) 566 618 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 619 END_2D … … 569 621 ! 570 622 DO jl = 1, jpl !-- first guess of tracer from v-flux 571 DO_2D( 0, 0, 0, 0)623 DO_2D( 0, 0, 1, 1 ) 572 624 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 573 625 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 576 628 END_2D 577 629 END DO 578 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )579 630 ! 580 631 DO jl = 1, jpl !-- flux in x-direction 581 DO_2D( 1, 0, 1, 0 )632 DO_2D( 0, 0, 1, 0 ) 582 633 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 634 END_2D … … 628 679 ! 629 680 DO jl = 1, jpl 630 DO_2D( 1, 0, 1, 0 )681 DO_2D( 1, 1, 1, 0 ) 631 682 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 683 END_2D 684 DO_2D( 1, 0, 1, 1 ) 632 685 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 633 686 END_2D … … 646 699 ! 647 700 DO jl = 1, jpl !-- flux in x-direction 648 DO_2D( 1, 0, 1, 0 )701 DO_2D( 1, 1, 1, 0 ) 649 702 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 650 703 END_2D … … 653 706 654 707 DO jl = 1, jpl !-- first guess of tracer from u-flux 655 DO_2D( 0, 0, 0, 0 )708 DO_2D( 1, 1, 0, 0 ) 656 709 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 657 710 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 660 713 END_2D 661 714 END DO 662 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )663 715 664 716 DO jl = 1, jpl !-- flux in y-direction 665 DO_2D( 1, 0, 1, 0 )717 DO_2D( 1, 0, 0, 0 ) 666 718 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 667 719 END_2D … … 672 724 ! 673 725 DO jl = 1, jpl !-- flux in y-direction 674 DO_2D( 1, 0, 1, 0)726 DO_2D( 1, 0, 1, 1 ) 675 727 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 676 728 END_2D … … 679 731 ! 680 732 DO jl = 1, jpl !-- first guess of tracer from v-flux 681 DO_2D( 0, 0, 0, 0)733 DO_2D( 0, 0, 1, 1 ) 682 734 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 683 735 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 686 738 END_2D 687 739 END DO 688 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )689 740 ! 690 741 DO jl = 1, jpl !-- flux in x-direction 691 DO_2D( 1, 0, 1, 0 )742 DO_2D( 0, 0, 1, 0 ) 692 743 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 693 744 END_2D … … 846 897 ! 847 898 DO jl = 1, jpl 848 DO_2D( 1, 0, 1, 0 )899 DO_2D( 0, 0, 1, 0 ) 849 900 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 850 901 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 855 906 ! 856 907 DO jl = 1, jpl 857 DO_2D( 1, 0, 1, 0 )908 DO_2D( 0, 0, 1, 0 ) 858 909 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 859 910 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 865 916 ! 866 917 DO jl = 1, jpl 867 DO_2D( 1, 0, 1, 0 )918 DO_2D( 0, 0, 1, 0 ) 868 919 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 869 920 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 879 930 ! 880 931 DO jl = 1, jpl 881 DO_2D( 1, 0, 1, 0 )932 DO_2D( 0, 0, 1, 0 ) 882 933 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 883 934 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 893 944 ! 894 945 DO jl = 1, jpl 895 DO_2D( 1, 0, 1, 0 )946 DO_2D( 0, 0, 1, 0 ) 896 947 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 897 948 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 914 965 IF( ll_neg ) THEN 915 966 DO jl = 1, jpl 916 DO_2D( 1, 0, 1, 0 )967 DO_2D( 0, 0, 1, 0 ) 917 968 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 918 969 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 924 975 ! !-- High order flux in i-direction --! 925 976 DO jl = 1, jpl 926 DO_2D( 1, 0, 1, 0 )977 DO_2D( 0, 0, 1, 0 ) 927 978 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 928 979 END_2D … … 957 1008 ! !-- Laplacian in j-direction --! 958 1009 DO jl = 1, jpl 959 DO_2D( 1, 0, 0, 0 ) 1010 DO_2D( 1, 0, 0, 0 ) ! First derivative (gradient) 960 1011 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 961 1012 END_2D 962 DO_2D( 0, 0, 0, 0 ) 1013 DO_2D( 0, 0, 0, 0 ) ! Second derivative (Laplacian) 963 1014 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 964 1015 END_2D … … 968 1019 ! !-- BiLaplacian in j-direction --! 969 1020 DO jl = 1, jpl 970 DO_2D( 1, 0, 0, 0 ) 1021 DO_2D( 1, 0, 0, 0 ) ! First derivative 971 1022 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 972 1023 END_2D 973 DO_2D( 0, 0, 0, 0 ) 1024 DO_2D( 0, 0, 0, 0 ) ! Second derivative 974 1025 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 975 1026 END_2D … … 982 1033 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 983 1034 DO jl = 1, jpl 984 DO_2D( 1, 0, 1, 0 )1035 DO_2D( 1, 0, 0, 0 ) 985 1036 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 986 1037 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 990 1041 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 991 1042 DO jl = 1, jpl 992 DO_2D( 1, 0, 1, 0 )1043 DO_2D( 1, 0, 0, 0 ) 993 1044 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 994 1045 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 999 1050 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1000 1051 DO jl = 1, jpl 1001 DO_2D( 1, 0, 1, 0 )1052 DO_2D( 1, 0, 0, 0 ) 1002 1053 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1003 1054 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1012 1063 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1013 1064 DO jl = 1, jpl 1014 DO_2D( 1, 0, 1, 0 )1065 DO_2D( 1, 0, 0, 0 ) 1015 1066 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1016 1067 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1025 1076 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1026 1077 DO jl = 1, jpl 1027 DO_2D( 1, 0, 1, 0 )1078 DO_2D( 1, 0, 0, 0 ) 1028 1079 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1029 1080 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1046 1097 IF( ll_neg ) THEN 1047 1098 DO jl = 1, jpl 1048 DO_2D( 1, 0, 1, 0 )1099 DO_2D( 1, 0, 0, 0 ) 1049 1100 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1050 1101 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1056 1107 ! !-- High order flux in j-direction --! 1057 1108 DO jl = 1, jpl 1058 DO_2D( 1, 0, 1, 0 )1109 DO_2D( 1, 0, 0, 0 ) 1059 1110 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1060 1111 END_2D … … 1092 1143 ! -------------------------------------------------- 1093 1144 DO jl = 1, jpl 1094 DO_2D( 1, 0, 1, 0 )1145 DO_2D( 0, 0, 1, 0 ) 1095 1146 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1147 END_2D 1148 DO_2D( 1, 0, 0, 0 ) 1096 1149 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1097 1150 END_2D … … 1199 1252 ! --------------------------------- 1200 1253 DO jl = 1, jpl 1201 DO_2D( 1, 0, 1, 0 )1254 DO_2D( 0, 0, 1, 0 ) 1202 1255 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1203 1256 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) … … 1210 1263 END_2D 1211 1264 1212 DO_2D( 1, 0, 1, 0 )1265 DO_2D( 1, 0, 0, 0 ) 1213 1266 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1214 1267 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) … … 1409 1462 1410 1463 1411 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 1464 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 1465 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 1412 1466 !!------------------------------------------------------------------- 1413 1467 !! *** ROUTINE Hbig *** … … 1423 1477 !! ** input : Max thickness of the surrounding 9-points 1424 1478 !!------------------------------------------------------------------- 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 1479 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1480 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max, psi_max ! max ice thick from surrounding 9-pts 1481 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pes_max 1482 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pei_max 1483 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 1428 1484 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 1485 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1486 ! 1487 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1488 REAL(wp) :: z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 1432 1489 !!------------------------------------------------------------------- 1433 1490 ! … … 1435 1492 ! 1436 1493 DO jl = 1, jpl 1437 1438 1494 DO_2D( 1, 1, 1, 1 ) 1439 1495 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN … … 1441 1497 ! ! -- check h_ip -- ! 1442 1498 ! 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 ) THEN1499 IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1444 1500 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1445 1501 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN … … 1468 1524 ENDIF 1469 1525 ! 1526 ! ! -- check s_i -- ! 1527 ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 1528 zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 1529 IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1530 zfra = psi_max(ji,jj,jl) / zsi 1531 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 1532 psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 1533 ENDIF 1534 ! 1470 1535 ENDIF 1471 1536 END_2D 1472 1537 END DO 1538 ! 1539 ! ! -- check e_i/v_i -- ! 1540 DO jl = 1, jpl 1541 DO_3D( 1, 1, 1, 1, 1, nlay_i ) 1542 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1543 ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 1544 zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 1545 IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1546 zfra = pei_max(ji,jj,jk,jl) / zei 1547 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1548 pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 1549 ENDIF 1550 ENDIF 1551 END_3D 1552 END DO 1553 ! ! -- check e_s/v_s -- ! 1554 DO jl = 1, jpl 1555 DO_3D( 1, 1, 1, 1, 1, nlay_s ) 1556 IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 1557 ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 1558 zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 1559 IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1560 zfra = pes_max(ji,jj,jk,jl) / zes 1561 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1562 pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 1563 ENDIF 1564 ENDIF 1565 END_3D 1566 END DO 1473 1567 ! 1474 1568 END SUBROUTINE Hbig … … 1526 1620 END SUBROUTINE Hsnow 1527 1621 1622 SUBROUTINE icemax3D( pice , pmax ) 1623 !!--------------------------------------------------------------------- 1624 !! *** ROUTINE icemax3D *** 1625 !! ** Purpose : compute the max of the 9 points around 1626 !!---------------------------------------------------------------------- 1627 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1628 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1629 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1630 INTEGER :: ji, jj, jl ! dummy loop indices 1631 !!---------------------------------------------------------------------- 1632 DO jl = 1, jpl 1633 DO jj = Njs0-1, Nje0+1 1634 DO ji = Nis0, Nie0 1635 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1636 END DO 1637 END DO 1638 DO jj = Njs0, Nje0 1639 DO ji = Nis0, Nie0 1640 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1641 END DO 1642 END DO 1643 END DO 1644 END SUBROUTINE icemax3D 1645 1646 SUBROUTINE icemax4D( pice , pmax ) 1647 !!--------------------------------------------------------------------- 1648 !! *** ROUTINE icemax4D *** 1649 !! ** Purpose : compute the max of the 9 points around 1650 !!---------------------------------------------------------------------- 1651 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1652 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1653 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1654 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1655 !!---------------------------------------------------------------------- 1656 jlay = SIZE( pice , 3 ) ! size of input arrays 1657 DO jl = 1, jpl 1658 DO jk = 1, jlay 1659 DO jj = Njs0-1, Nje0+1 1660 DO ji = Nis0, Nie0 1661 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1662 END DO 1663 END DO 1664 DO jj = Njs0, Nje0 1665 DO ji = Nis0, Nie0 1666 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1667 END DO 1668 END DO 1669 END DO 1670 END DO 1671 END SUBROUTINE icemax4D 1528 1672 1529 1673 #else
Note: See TracChangeset
for help on using the changeset viewer.
