- Timestamp:
- 2020-11-30T12:29:45+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion
- Files:
-
- 17 edited
-
. (modified) (1 prop)
-
src/ICE/ice.F90 (modified) (2 diffs)
-
src/ICE/ice1d.F90 (modified) (2 diffs)
-
src/ICE/icecor.F90 (modified) (3 diffs)
-
src/ICE/icectl.F90 (modified) (7 diffs)
-
src/ICE/icedyn_adv_pra.F90 (modified) (19 diffs)
-
src/ICE/icedyn_adv_umx.F90 (modified) (44 diffs)
-
src/ICE/icedyn_rdgrft.F90 (modified) (1 diff)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (11 diffs)
-
src/ICE/iceitd.F90 (modified) (2 diffs)
-
src/ICE/icestp.F90 (modified) (11 diffs)
-
src/ICE/icetab.F90 (modified) (4 diffs)
-
src/ICE/icethd.F90 (modified) (13 diffs)
-
src/ICE/icethd_dh.F90 (modified) (2 diffs)
-
src/ICE/icethd_do.F90 (modified) (4 diffs)
-
src/ICE/icethd_ent.F90 (modified) (2 diffs)
-
src/ICE/iceupdate.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette@135 07sette10 ^/utils/CI/sette@13559 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/ice.F90
r13472 r13917 391 391 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 392 392 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_aice !: ice conc. variation [s-1] 394 ! 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_mass !: advection of mass (kg/m2/s) 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_salt !: advection of salt (g/m2/s) 397 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_heat !: advection of heat (W/m2) 393 398 ! 394 399 !!---------------------------------------------------------------------- … … 493 498 ! * Ice diagnostics 494 499 ii = ii + 1 495 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 496 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 497 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), STAT=ierr(ii) ) 500 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 501 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 502 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), & 503 & diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 498 504 499 505 ! * Ice conservation -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/ice1d.F90
r13472 r13917 145 145 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sst_1d 146 146 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sss_1d 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: frq_m_1d 147 148 148 149 ! convergence check … … 225 226 ! 226 227 ii = ii + 1 227 ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) )228 ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , frq_m_1d(jpij) , STAT=ierr(ii) ) 228 229 ! 229 230 ii = ii + 1 -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icecor.F90
r13497 r13917 55 55 INTEGER :: ji, jj, jk, jl ! dummy loop indices 56 56 REAL(wp) :: zsal, zzc 57 REAL(wp), DIMENSION(jpi,jpj) :: zafx ! concentration trends diag58 57 !!---------------------------------------------------------------------- 59 58 ! controls … … 95 94 zsal = sv_i(ji,jj,jl) 96 95 sv_i(ji,jj,jl) = MIN( MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl) ) 97 sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 96 IF( kn /= 0 ) & ! no ice-ocean exchanges if kn=0 (for bdy for instance) otherwise conservation diags will fail 97 & sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 98 98 END_2D 99 99 END DO 100 100 ENDIF 101 ! !-----------------------------------------------------102 CALL ice_var_zapsmall ! Zap small values !103 ! !-----------------------------------------------------104 101 102 IF( kn /= 0 ) THEN ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx) 103 ! otherwise conservation diags will fail 104 ! !----------------------------------------------------- 105 CALL ice_var_zapsmall ! Zap small values ! 106 ! !----------------------------------------------------- 107 ENDIF 105 108 ! !----------------------------------------------------- 106 109 IF( kn == 2 ) THEN ! Ice drift case: Corrections to avoid wrong values ! … … 115 118 CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 116 119 ENDIF 117 118 ! !-----------------------------------------------------119 SELECT CASE( kn ) ! Diagnostics !120 ! !-----------------------------------------------------121 CASE( 1 ) !--- dyn trend diagnostics122 !123 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN124 diag_heat(:,:) = - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & ! W.m-2125 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice126 diag_sice(:,:) = SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi127 diag_vice(:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi128 diag_vsnw(:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos129 ENDIF130 ! ! concentration tendency (dynamics)131 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN132 zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice133 CALL iom_put( 'afxdyn' , zafx )134 ENDIF135 !136 CASE( 2 ) !--- thermo trend diagnostics & ice aging137 !138 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! ice natural aging incrementation139 !140 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN141 diag_heat(:,:) = diag_heat(:,:) &142 & - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &143 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice144 diag_sice(:,:) = diag_sice(:,:) &145 & + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi146 diag_vice(:,:) = diag_vice(:,:) &147 & + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi148 diag_vsnw(:,:) = diag_vsnw(:,:) &149 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos150 CALL iom_put ( 'hfxdhc' , diag_heat )151 ENDIF152 ! ! concentration tendency (total + thermo)153 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN154 zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice155 CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice )156 CALL iom_put( 'afxtot' , zafx )157 ENDIF158 !159 END SELECT160 120 ! 161 121 ! controls -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icectl.F90
r13472 r13917 43 43 PUBLIC ice_prt 44 44 PUBLIC ice_prt3D 45 PUBLIC ice_drift_wri 46 PUBLIC ice_drift_init 45 47 46 48 ! thresold rates for conservation … … 49 51 REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 50 52 REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 53 54 ! for drift outputs 55 CHARACTER(LEN=50) :: clname="icedrift_diagnostics.ascii" ! ascii filename 56 INTEGER :: numicedrift ! outfile unit 57 REAL(wp) :: rdiag_icemass, rdiag_icesalt, rdiag_iceheat 58 REAL(wp) :: rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat 51 59 52 60 !! * Substitutions … … 132 140 133 141 ! -- advection scheme is conservative? -- ! 134 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather)135 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) ! must be close to 0 (only for Prather)142 zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 143 zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 136 144 137 145 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 156 164 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_amax 157 165 ! check if advection scheme is conservative 158 ! only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 159 ! so the formulation for conservation is different (and not coded) 160 ! it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 161 !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 162 ! & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 166 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 167 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 168 IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 169 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * rdt_ice 163 170 ENDIF 164 171 ! … … 186 193 ! water flux 187 194 ! -- mass diag -- ! 188 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 195 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 196 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 189 197 190 198 ! -- salt diag -- ! 191 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t )199 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 192 200 193 201 ! -- heat diag -- ! 194 ! clem: not the good formulation 195 !!zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr & 196 !! & ) * e1e2t ) 202 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 203 ! equivalent to this: 204 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 205 !! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 206 !! & ) * e1e2t ) 197 207 198 208 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 204 214 IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 205 215 & WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zdiag_salt * rDt_ice 206 !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 216 IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 217 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 207 218 ENDIF 208 219 ! … … 725 736 726 737 END SUBROUTINE ice_prt3D 738 739 740 SUBROUTINE ice_drift_wri( kt ) 741 !!------------------------------------------------------------------- 742 !! *** ROUTINE ice_drift_wri *** 743 !! 744 !! ** Purpose : conservation of mass, salt and heat 745 !! write the drift in a ascii file at each time step 746 !! and the total run drifts 747 !!------------------------------------------------------------------- 748 INTEGER, INTENT(in) :: kt ! ice time-step index 749 ! 750 INTEGER :: ji, jj 751 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 752 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 753 !!------------------------------------------------------------------- 754 ! 755 IF( kt == nit000 .AND. lwp ) THEN 756 WRITE(numout,*) 757 WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 758 WRITE(numout,*) '~~~~~~~~~~~~~' 759 ENDIF 760 ! 761 ! 2D budgets (must be close to 0) 762 IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 763 DO_2D( 1, 1, 1, 1 ) 764 zdiag_mass2D(ji,jj) = wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 765 & + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 766 zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 767 zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 768 END_2D 769 ! 770 ! write outputs 771 CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 772 CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 773 CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 774 ENDIF 775 776 ! -- mass diag -- ! 777 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 778 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 779 zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 780 781 ! -- salt diag -- ! 782 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 783 zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 784 785 ! -- heat diag -- ! 786 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 787 zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 788 789 ! ! write out to file 790 IF( lwp ) THEN 791 ! check global drift (must be close to 0) 792 WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift [kg]', zdiag_mass 793 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift [kg]', zdiag_salt 794 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift [W] ', zdiag_heat 795 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 796 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'mass drift adv [kg]', zdiag_adv_mass 797 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift adv [kg]', zdiag_adv_salt 798 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift adv [W] ', zdiag_adv_heat 799 ENDIF 800 ! ! drifts 801 rdiag_icemass = rdiag_icemass + zdiag_mass 802 rdiag_icesalt = rdiag_icesalt + zdiag_salt 803 rdiag_iceheat = rdiag_iceheat + zdiag_heat 804 rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 805 rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 806 rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 807 ! 808 ! ! output drifts and close ascii file 809 IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 810 ! to ascii file 811 WRITE(numicedrift,*) '******************************************' 812 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift [kg]', rdiag_icemass 813 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 814 WRITE(numicedrift,*) '******************************************' 815 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift [kg]', rdiag_icesalt 816 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 817 WRITE(numicedrift,*) '******************************************' 818 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift [W] ', rdiag_iceheat 819 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 820 CLOSE( numicedrift ) 821 ! 822 ! to ocean output 823 WRITE(numout,*) 824 WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 825 WRITE(numout,*) '~~~~~~~~~~~~~' 826 ! check global drift (must be close to 0) 827 WRITE(numout,*) ' sea-ice mass drift [kg] = ', rdiag_icemass 828 WRITE(numout,*) ' sea-ice salt drift [kg] = ', rdiag_icesalt 829 WRITE(numout,*) ' sea-ice heat drift [W] = ', rdiag_iceheat 830 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 831 WRITE(numout,*) ' sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 832 WRITE(numout,*) ' sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 833 WRITE(numout,*) ' sea-ice heat drift adv [W] = ', rdiag_adv_iceheat 834 ENDIF 835 ! 836 END SUBROUTINE ice_drift_wri 837 838 SUBROUTINE ice_drift_init 839 !!---------------------------------------------------------------------- 840 !! *** ROUTINE ice_drift_init *** 841 !! 842 !! ** Purpose : create output file, initialise arrays 843 !!---------------------------------------------------------------------- 844 ! 845 IF( .NOT.ln_icediachk ) RETURN ! exit 846 ! 847 IF(lwp) THEN 848 WRITE(numout,*) 849 WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 850 WRITE(numout,*) '~~~~~~~~~~~~~' 851 WRITE(numout,*) 852 ! 853 ! create output ascii file 854 CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 855 WRITE(numicedrift,*) 'Timestep Drifts' 856 WRITE(numicedrift,*) '******************************************' 857 ENDIF 858 ! 859 rdiag_icemass = 0._wp 860 rdiag_icesalt = 0._wp 861 rdiag_iceheat = 0._wp 862 rdiag_adv_icemass = 0._wp 863 rdiag_adv_icesalt = 0._wp 864 rdiag_adv_iceheat = 0._wp 865 ! 866 END SUBROUTINE ice_drift_init 727 867 728 868 #else -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icedyn_adv_pra.F90
r13497 r13917 88 88 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 89 89 INTEGER :: icycle ! number of sub-timestep for the advection 90 REAL(wp) :: zdt 90 REAL(wp) :: zdt, z1_dt ! - - 91 91 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 92 92 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 … … 100 100 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 101 101 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei 102 !! diagnostics 103 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 102 104 !!---------------------------------------------------------------------- 103 105 ! … … 109 111 ELSEWHERE ; zs_i(:,:,:) = 0._wp 110 112 END WHERE 111 DO jl = 1, jpl 112 DO_2D( 0, 0, 0, 0 ) 113 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 114 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 115 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 116 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 117 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 118 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 119 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 120 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 121 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 122 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 123 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 124 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 125 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 126 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 127 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 128 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 129 END_2D 130 END DO 113 CALL icemax3D( ph_i , zhi_max ) 114 CALL icemax3D( ph_s , zhs_max ) 115 CALL icemax3D( ph_ip, zhip_max) 116 CALL icemax3D( zs_i , zsi_max ) 131 117 CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 132 118 ! … … 141 127 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 142 128 END WHERE 143 END DO 144 DO jl = 1, jpl 145 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 146 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 147 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 148 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 149 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 150 END_3D 151 END DO 152 DO jl = 1, jpl 153 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 154 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 155 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 156 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 157 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 158 END_3D 159 END DO 160 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 161 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 129 END DO 130 CALL icemax4D( ze_i , zei_max ) 131 CALL icemax4D( ze_s , zes_max ) 132 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 133 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 162 134 ! 163 135 ! … … 175 147 ENDIF 176 148 zdt = rDt_ice / REAL(icycle) 149 z1_dt = 1._wp / zdt 177 150 178 151 ! --- transport --- ! … … 181 154 182 155 DO jt = 1, icycle 156 157 ! diagnostics 158 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 159 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 160 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 161 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 183 162 184 163 ! record at_i before advection (for open water) … … 279 258 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 280 259 ENDIF 281 ENDIF 282 ! 260 ENDIF 261 ! 262 ENDIF 263 264 ! --- Lateral boundary conditions --- ! 265 ! caution: for gradients (sx and sy) the sign changes 266 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp & ! ice volume 267 & , sxxice, 'T', 1._wp, syyice, 'T', 1._wp, sxyice, 'T', 1._wp & 268 & , z0snw , 'T', 1._wp, sxsn , 'T', -1._wp, sysn , 'T', -1._wp & ! snw volume 269 & , sxxsn , 'T', 1._wp, syysn , 'T', 1._wp, sxysn , 'T', 1._wp ) 270 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp & ! ice salinity 271 & , sxxsal, 'T', 1._wp, syysal, 'T', 1._wp, sxysal, 'T', 1._wp & 272 & , z0ai , 'T', 1._wp, sxa , 'T', -1._wp, sya , 'T', -1._wp & ! ice concentration 273 & , sxxa , 'T', 1._wp, syya , 'T', 1._wp, sxya , 'T', 1._wp ) 274 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age 275 & , sxxage, 'T', 1._wp, syyage, 'T', 1._wp, sxyage, 'T', 1._wp ) 276 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es , 'T', 1._wp, sxc0 , 'T', -1._wp, syc0 , 'T', -1._wp & ! snw enthalpy 277 & , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp ) 278 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy 279 & , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp ) 280 IF ( ln_pnd_LEV ) THEN 281 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction 282 & , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp & 283 & , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume 284 & , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp ) 285 IF ( ln_pnd_lids ) THEN 286 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp & ! melt pond lid volume 287 & , sxxvl,'T', 1._wp, syyvl,'T', 1._wp, sxyvl,'T', 1._wp ) 288 ENDIF 283 289 ENDIF 284 290 … … 312 318 END_2D 313 319 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp ) 320 ! 321 ! --- diagnostics --- ! 322 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 323 & - zdiag_adv_mass(:,:) ) * z1_dt 324 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 325 & - zdiag_adv_salt(:,:) ) * z1_dt 326 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 327 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 328 & - zdiag_adv_heat(:,:) ) * z1_dt 314 329 ! 315 330 ! --- Ensure non-negative fields --- ! … … 350 365 !! 351 366 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 367 INTEGER :: jj0 ! dummy loop indices 352 368 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 353 369 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 354 370 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 371 REAL(wp) :: zpsm, zps0 372 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 355 373 REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace 356 374 REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - - 357 375 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 358 376 !----------------------------------------------------------------------- 377 ! in order to avoid lbc_lnk (communications): 378 ! jj loop must be 1:jpj if adv_x is called first 379 ! and 2:jpj-1 if adv_x is called second 380 jj0 = NINT(pcrh) 359 381 ! 360 382 jcat = SIZE( ps0 , 3 ) ! size of input arrays … … 363 385 ! 364 386 ! Limitation of moments. 365 DO_2D( 0, 0, 1, 1 ) 366 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 367 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 368 ! 369 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 370 zs1max = 1.5 * zslpmax 371 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 372 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 373 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 374 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 375 376 ps0 (ji,jj,jl) = zslpmax 377 psx (ji,jj,jl) = zs1new * rswitch 378 psxx(ji,jj,jl) = zs2new * rswitch 379 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 380 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 381 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 382 END_2D 383 384 ! Calculate fluxes and moments between boxes i<-->i+1 385 DO_2D( 0, 0, 1, 1 ) ! Flux from i to i+1 WHEN u GT 0 386 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 387 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 388 zalfq = zalf * zalf 389 zalf1 = 1.0 - zalf 390 zalf1q = zalf1 * zalf1 391 ! 392 zfm (ji,jj) = zalf * psm (ji,jj,jl) 393 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 394 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 395 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 396 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 397 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 398 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 399 400 ! Readjust moments remaining in the box. 401 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 402 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 403 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 404 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 405 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 406 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 407 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 408 END_2D 409 410 DO_2D( 0, 0, 1, 0 ) ! Flux from i+1 to i when u LT 0. 411 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 412 zalg (ji,jj) = zalf 413 zalfq = zalf * zalf 414 zalf1 = 1.0 - zalf 415 zalg1 (ji,jj) = zalf1 416 zalf1q = zalf1 * zalf1 417 zalg1q(ji,jj) = zalf1q 418 ! 419 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 420 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 421 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 422 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 423 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 424 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 425 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 426 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 427 END_2D 428 429 DO_2D( 0, 0, 0, 0 ) ! Readjust moments remaining in the box. 430 zbt = zbet(ji-1,jj) 431 zbt1 = 1.0 - zbet(ji-1,jj) 432 ! 433 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 434 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 435 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 436 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 437 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 438 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 439 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 440 END_2D 441 442 ! Put the temporary moments into appropriate neighboring boxes. 443 DO_2D( 0, 0, 0, 0 ) ! Flux from i to i+1 IF u GT 0. 444 zbt = zbet(ji-1,jj) 445 zbt1 = 1.0 - zbet(ji-1,jj) 446 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 447 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 448 zalf1 = 1.0 - zalf 449 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 450 ! 451 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 452 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 453 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 454 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 455 & + zbt1 * psxx(ji,jj,jl) 456 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 457 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 458 & + zbt1 * psxy(ji,jj,jl) 459 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 460 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 461 END_2D 462 463 DO_2D( 0, 0, 0, 0 ) ! Flux from i+1 to i IF u LT 0. 464 zbt = zbet(ji,jj) 465 zbt1 = 1.0 - zbet(ji,jj) 466 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 467 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 468 zalf1 = 1.0 - zalf 469 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 470 ! 471 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 472 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 473 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 474 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 475 & + ( zalf1 - zalf ) * ztemp ) ) 476 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 477 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 478 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 479 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 480 END_2D 481 387 DO jj = Njs0 - jj0, Nje0 + jj0 388 389 DO ji = Nis0 - 1, Nie0 + 1 390 391 zpsm = psm (ji,jj,jl) ! optimization 392 zps0 = ps0 (ji,jj,jl) 393 zpsx = psx (ji,jj,jl) 394 zpsxx = psxx(ji,jj,jl) 395 zpsy = psy (ji,jj,jl) 396 zpsyy = psyy(ji,jj,jl) 397 zpsxy = psxy(ji,jj,jl) 398 399 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 400 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 401 ! 402 zslpmax = MAX( 0._wp, zps0 ) 403 zs1max = 1.5 * zslpmax 404 zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) ) 405 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 406 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 407 408 zps0 = zslpmax 409 zpsx = zs1new * rswitch 410 zpsxx = zs2new * rswitch 411 zpsy = zpsy * rswitch 412 zpsyy = zpsyy * rswitch 413 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 414 415 ! Calculate fluxes and moments between boxes i<-->i+1 416 ! ! Flux from i to i+1 WHEN u GT 0 417 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 418 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 419 zalfq = zalf * zalf 420 zalf1 = 1.0 - zalf 421 zalf1q = zalf1 * zalf1 422 ! 423 zfm (ji,jj) = zalf * zpsm 424 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 425 zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx ) 426 zfxx(ji,jj) = zalf * zpsxx * zalfq 427 zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy ) 428 zfxy(ji,jj) = zalfq * zpsxy 429 zfyy(ji,jj) = zalf * zpsyy 430 431 ! ! Readjust moments remaining in the box. 432 zpsm = zpsm - zfm(ji,jj) 433 zps0 = zps0 - zf0(ji,jj) 434 zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 435 zpsxx = zalf1 * zalf1q * zpsxx 436 zpsy = zpsy - zfy (ji,jj) 437 zpsyy = zpsyy - zfyy(ji,jj) 438 zpsxy = zalf1q * zpsxy 439 ! 440 psm (ji,jj,jl) = zpsm ! optimization 441 ps0 (ji,jj,jl) = zps0 442 psx (ji,jj,jl) = zpsx 443 psxx(ji,jj,jl) = zpsxx 444 psy (ji,jj,jl) = zpsy 445 psyy(ji,jj,jl) = zpsyy 446 psxy(ji,jj,jl) = zpsxy 447 ! 448 END DO 449 450 DO ji = Nis0 - 1, Nie0 451 ! ! Flux from i+1 to i when u LT 0. 452 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 453 zalg (ji,jj) = zalf 454 zalfq = zalf * zalf 455 zalf1 = 1.0 - zalf 456 zalg1 (ji,jj) = zalf1 457 zalf1q = zalf1 * zalf1 458 zalg1q(ji,jj) = zalf1q 459 ! 460 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 461 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 462 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 463 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 464 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 465 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 466 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 467 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 468 END DO 469 470 DO ji = Nis0, Nie0 471 ! 472 zpsm = psm (ji,jj,jl) ! optimization 473 zps0 = ps0 (ji,jj,jl) 474 zpsx = psx (ji,jj,jl) 475 zpsxx = psxx(ji,jj,jl) 476 zpsy = psy (ji,jj,jl) 477 zpsyy = psyy(ji,jj,jl) 478 zpsxy = psxy(ji,jj,jl) 479 ! ! Readjust moments remaining in the box. 480 zbt = zbet(ji-1,jj) 481 zbt1 = 1.0 - zbet(ji-1,jj) 482 ! 483 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 484 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 485 zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 486 zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 487 zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) ) 488 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 489 zpsxy = zalg1q(ji-1,jj) * zpsxy 490 491 ! Put the temporary moments into appropriate neighboring boxes. 492 ! ! Flux from i to i+1 IF u GT 0. 493 zbt = zbet(ji-1,jj) 494 zbt1 = 1.0 - zbet(ji-1,jj) 495 zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 496 zalf = zbt * zfm(ji-1,jj) / zpsm 497 zalf1 = 1.0 - zalf 498 ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 499 ! 500 zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 501 zpsx = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 502 zpsxx = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx & 503 & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 504 & + zbt1 * zpsxx 505 zpsxy = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy & 506 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * zpsy ) ) & 507 & + zbt1 * zpsxy 508 zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy 509 zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 510 511 ! ! Flux from i+1 to i IF u LT 0. 512 zbt = zbet(ji,jj) 513 zbt1 = 1.0 - zbet(ji,jj) 514 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 515 zalf = zbt1 * zfm(ji,jj) / zpsm 516 zalf1 = 1.0 - zalf 517 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 518 ! 519 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 520 zpsx = zbt * zpsx + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 521 zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 522 & + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) & 523 & + ( zalf1 - zalf ) * ztemp ) ) 524 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 525 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 526 zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) ) 527 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 528 ! 529 psm (ji,jj,jl) = zpsm ! optimization 530 ps0 (ji,jj,jl) = zps0 531 psx (ji,jj,jl) = zpsx 532 psxx(ji,jj,jl) = zpsxx 533 psy (ji,jj,jl) = zpsy 534 psyy(ji,jj,jl) = zpsyy 535 psxy(ji,jj,jl) = zpsxy 536 END DO 537 ! 538 END DO 539 ! 482 540 END DO 483 484 !-- Lateral boundary conditions 485 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 486 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 487 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 488 ! 541 ! 489 542 END SUBROUTINE adv_x 490 543 … … 507 560 !! 508 561 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 562 INTEGER :: ji0 ! dummy loop indices 509 563 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 510 564 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 511 565 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 566 REAL(wp) :: zpsm, zps0 567 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 512 568 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 513 569 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - 514 570 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 515 571 !--------------------------------------------------------------------- 572 ! in order to avoid lbc_lnk (communications): 573 ! ji loop must be 1:jpi if adv_y is called first 574 ! and 2:jpi-1 if adv_y is called second 575 ji0 = NINT(pcrh) 516 576 ! 517 577 jcat = SIZE( ps0 , 3 ) ! size of input arrays … … 520 580 ! 521 581 ! Limitation of moments. 522 DO_2D( 1, 1, 0, 0 ) 523 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 524 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 525 ! 526 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 582 DO_2D( 1, 1, ji0, ji0 ) 583 ! 584 zpsm = psm (ji,jj,jl) ! optimization 585 zps0 = ps0 (ji,jj,jl) 586 zpsx = psx (ji,jj,jl) 587 zpsxx = psxx(ji,jj,jl) 588 zpsy = psy (ji,jj,jl) 589 zpsyy = psyy(ji,jj,jl) 590 zpsxy = psxy(ji,jj,jl) 591 ! 592 ! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 593 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 594 ! 595 zslpmax = MAX( 0._wp, zps0 ) 527 596 zs1max = 1.5 * zslpmax 528 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 529 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 530 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 597 zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) ) 598 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 531 599 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 532 600 ! 533 ps0 (ji,jj,jl) = zslpmax 534 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 535 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 536 psy (ji,jj,jl) = zs1new * rswitch 537 psyy(ji,jj,jl) = zs2new * rswitch 538 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 539 END_2D 540 541 ! Calculate fluxes and moments between boxes j<-->j+1 542 DO_2D( 1, 1, 0, 0 ) ! Flux from j to j+1 WHEN v GT 0 601 zps0 = zslpmax 602 zpsx = zpsx * rswitch 603 zpsxx = zpsxx * rswitch 604 zpsy = zs1new * rswitch 605 zpsyy = zs2new * rswitch 606 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 607 608 ! Calculate fluxes and moments between boxes j<-->j+1 609 ! ! Flux from j to j+1 WHEN v GT 0 543 610 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 544 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)611 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 545 612 zalfq = zalf * zalf 546 613 zalf1 = 1.0 - zalf 547 614 zalf1q = zalf1 * zalf1 548 615 ! 549 zfm (ji,jj) = zalf * psm(ji,jj,jl) 550 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 551 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 552 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 553 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 554 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 555 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 556 ! 557 ! Readjust moments remaining in the box. 558 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 559 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 560 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 561 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 562 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 563 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 564 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 616 zfm (ji,jj) = zalf * zpsm 617 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) ) 618 zfy (ji,jj) = zalfq *( zpsy + 3.0*zalf1*zpsyy ) 619 zfyy(ji,jj) = zalf * zalfq * zpsyy 620 zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy ) 621 zfxy(ji,jj) = zalfq * zpsxy 622 zfxx(ji,jj) = zalf * zpsxx 623 ! 624 ! ! Readjust moments remaining in the box. 625 zpsm = zpsm - zfm(ji,jj) 626 zps0 = zps0 - zf0(ji,jj) 627 zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 628 zpsyy = zalf1 * zalf1q * zpsyy 629 zpsx = zpsx - zfx(ji,jj) 630 zpsxx = zpsxx - zfxx(ji,jj) 631 zpsxy = zalf1q * zpsxy 632 ! 633 psm (ji,jj,jl) = zpsm ! optimization 634 ps0 (ji,jj,jl) = zps0 635 psx (ji,jj,jl) = zpsx 636 psxx(ji,jj,jl) = zpsxx 637 psy (ji,jj,jl) = zpsy 638 psyy(ji,jj,jl) = zpsyy 639 psxy(ji,jj,jl) = zpsxy 565 640 END_2D 566 641 ! 567 DO_2D( 1, 0, 0, 0 ) ! Flux from j+1 to j when v LT 0. 642 DO_2D( 1, 0, ji0, ji0 ) 643 ! ! Flux from j+1 to j when v LT 0. 568 644 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 569 645 zalg (ji,jj) = zalf … … 584 660 END_2D 585 661 586 ! Readjust moments remaining in the box.587 DO_2D( 0, 0, 0, 0 )662 DO_2D( 0, 0, ji0, ji0 ) 663 ! ! Readjust moments remaining in the box. 588 664 zbt = zbet(ji,jj-1) 589 665 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 590 666 ! 591 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 592 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 593 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 594 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 595 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 596 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 597 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 667 zpsm = psm (ji,jj,jl) ! optimization 668 zps0 = ps0 (ji,jj,jl) 669 zpsx = psx (ji,jj,jl) 670 zpsxx = psxx(ji,jj,jl) 671 zpsy = psy (ji,jj,jl) 672 zpsyy = psyy(ji,jj,jl) 673 zpsxy = psxy(ji,jj,jl) 674 ! 675 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 676 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 677 zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 678 zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 679 zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) ) 680 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 681 zpsxy = zalg1q(ji,jj-1) * zpsxy 682 683 ! Put the temporary moments into appropriate neighboring boxes. 684 ! ! Flux from j to j+1 IF v GT 0. 685 zbt = zbet(ji,jj-1) 686 zbt1 = 1.0 - zbet(ji,jj-1) 687 zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm 688 zalf = zbt * zfm(ji,jj-1) / zpsm 689 zalf1 = 1.0 - zalf 690 ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 691 ! 692 zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 693 zpsy = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp ) & 694 & + zbt1 * zpsy 695 zpsyy = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy & 696 & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 697 & + zbt1 * zpsyy 698 zpsxy = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy & 699 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) & 700 & + zbt1 * zpsxy 701 zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx 702 zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 703 704 ! ! Flux from j+1 to j IF v LT 0. 705 zbt = zbet(ji,jj) 706 zbt1 = 1.0 - zbet(ji,jj) 707 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 708 zalf = zbt1 * zfm(ji,jj) / zpsm 709 zalf1 = 1.0 - zalf 710 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 711 ! 712 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 713 zpsy = zbt * zpsy + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 714 zpsyy = zbt * zpsyy + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 715 & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) & 716 & + ( zalf1 - zalf ) * ztemp ) ) 717 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 718 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 719 zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) ) 720 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 721 ! 722 psm (ji,jj,jl) = zpsm ! optimization 723 ps0 (ji,jj,jl) = zps0 724 psx (ji,jj,jl) = zpsx 725 psxx(ji,jj,jl) = zpsxx 726 psy (ji,jj,jl) = zpsy 727 psyy(ji,jj,jl) = zpsyy 728 psxy(ji,jj,jl) = zpsxy 598 729 END_2D 599 600 ! Put the temporary moments into appropriate neighboring boxes. 601 DO_2D( 0, 0, 0, 0 ) ! Flux from j to j+1 IF v GT 0. 602 zbt = zbet(ji,jj-1) 603 zbt1 = 1.0 - zbet(ji,jj-1) 604 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 605 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 606 zalf1 = 1.0 - zalf 607 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 608 ! 609 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 610 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 611 & + zbt1 * psy(ji,jj,jl) 612 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 613 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 614 & + zbt1 * psyy(ji,jj,jl) 615 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 616 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 617 & + zbt1 * psxy(ji,jj,jl) 618 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 619 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 620 END_2D 621 622 DO_2D( 0, 0, 0, 0 ) ! Flux from j+1 to j IF v LT 0. 623 zbt = zbet(ji,jj) 624 zbt1 = 1.0 - zbet(ji,jj) 625 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 626 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 627 zalf1 = 1.0 - zalf 628 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 629 ! 630 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 631 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 632 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 633 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 634 & + ( zalf1 - zalf ) * ztemp ) ) 635 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 636 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 637 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 638 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 639 END_2D 640 730 ! 641 731 END DO 642 643 !-- Lateral boundary conditions644 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp &645 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes646 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp )647 732 ! 648 733 END SUBROUTINE adv_y … … 872 957 ! 873 958 ! ! ice thickness 874 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice )875 CALL iom_get( numrir, jpdom_auto, 'syice' , syice )959 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp ) 960 CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp ) 876 961 CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 877 962 CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 878 963 CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 879 964 ! ! snow thickness 880 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn )881 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn )965 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn , psgn = -1._wp ) 966 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn , psgn = -1._wp ) 882 967 CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn ) 883 968 CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn ) 884 969 CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn ) 885 970 ! ! ice concentration 886 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa )887 CALL iom_get( numrir, jpdom_auto, 'sya' , sya )971 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa , psgn = -1._wp ) 972 CALL iom_get( numrir, jpdom_auto, 'sya' , sya , psgn = -1._wp ) 888 973 CALL iom_get( numrir, jpdom_auto, 'sxxa' , sxxa ) 889 974 CALL iom_get( numrir, jpdom_auto, 'syya' , syya ) 890 975 CALL iom_get( numrir, jpdom_auto, 'sxya' , sxya ) 891 976 ! ! ice salinity 892 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal )893 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal )977 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp ) 978 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp ) 894 979 CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 895 980 CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 896 981 CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 897 982 ! ! ice age 898 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage )899 CALL iom_get( numrir, jpdom_auto, 'syage' , syage )983 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp ) 984 CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp ) 900 985 CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 901 986 CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) … … 904 989 DO jk = 1, nlay_s 905 990 WRITE(zchar1,'(I2.2)') jk 906 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)907 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:)991 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxc0 (:,:,jk,:) = z3d(:,:,:) 992 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; syc0 (:,:,jk,:) = z3d(:,:,:) 908 993 znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:) 909 994 znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:) … … 913 998 DO jk = 1, nlay_i 914 999 WRITE(zchar1,'(I2.2)') jk 915 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:)916 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:)1000 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxe (:,:,jk,:) = z3d(:,:,:) 1001 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sye (:,:,jk,:) = z3d(:,:,:) 917 1002 znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:) 918 1003 znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:) … … 921 1006 ! 922 1007 IF( ln_pnd_LEV ) THEN ! melt pond fraction 923 IF( iom_varid( numr or, 'sxap', ldstop = .FALSE. ) > 0 ) THEN924 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap )925 CALL iom_get( numrir, jpdom_auto, 'syap' , syap )1008 IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 1009 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 1010 CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp ) 926 1011 CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 927 1012 CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 928 1013 CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 929 1014 ! ! melt pond volume 930 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp )931 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp )1015 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp ) 1016 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp ) 932 1017 CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 933 1018 CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) … … 939 1024 ! 940 1025 IF ( ln_pnd_lids ) THEN ! melt pond lid volume 941 IF( iom_varid( numr or, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN942 CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl )943 CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl )1026 IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 1027 CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp ) 1028 CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp ) 944 1029 CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 945 1030 CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) … … 1056 1141 END SUBROUTINE adv_pra_rst 1057 1142 1143 SUBROUTINE icemax3D( pice , pmax ) 1144 !!--------------------------------------------------------------------- 1145 !! *** ROUTINE icemax3D *** 1146 !! ** Purpose : compute the max of the 9 points around 1147 !!---------------------------------------------------------------------- 1148 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1149 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1150 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1151 INTEGER :: ji, jj, jl ! dummy loop indices 1152 !!---------------------------------------------------------------------- 1153 DO jl = 1, jpl 1154 DO jj = Njs0-1, Nje0+1 1155 DO ji = Nis0, Nie0 1156 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1157 END DO 1158 END DO 1159 DO jj = Njs0, Nje0 1160 DO ji = Nis0, Nie0 1161 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1162 END DO 1163 END DO 1164 END DO 1165 END SUBROUTINE icemax3D 1166 1167 SUBROUTINE icemax4D( pice , pmax ) 1168 !!--------------------------------------------------------------------- 1169 !! *** ROUTINE icemax4D *** 1170 !! ** Purpose : compute the max of the 9 points around 1171 !!---------------------------------------------------------------------- 1172 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1173 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1174 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1175 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1176 !!---------------------------------------------------------------------- 1177 jlay = SIZE( pice , 3 ) ! size of input arrays 1178 DO jl = 1, jpl 1179 DO jk = 1, jlay 1180 DO jj = Njs0-1, Nje0+1 1181 DO ji = Nis0, Nie0 1182 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1183 END DO 1184 END DO 1185 DO jj = Njs0, Nje0 1186 DO ji = Nis0, Nie0 1187 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1188 END DO 1189 END DO 1190 END DO 1191 END DO 1192 END SUBROUTINE icemax4D 1193 1058 1194 #else 1059 1195 !!---------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icedyn_adv_umx.F90
r13497 r13917 92 92 INTEGER :: icycle ! number of sub-timestep for the advection 93 93 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 94 REAL(wp) :: zdt, z vi_cen94 REAL(wp) :: zdt, z1_dt, zvi_cen 95 95 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 96 96 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box … … 104 104 ! 105 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 106 108 !!---------------------------------------------------------------------- 107 109 ! … … 113 115 ELSEWHERE ; zs_i(:,:,:) = 0._wp 114 116 END WHERE 115 DO jl = 1, jpl 116 DO_2D( 0, 0, 0, 0 ) 117 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 118 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 119 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 120 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 121 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 122 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 123 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 124 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 125 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 126 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 127 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 128 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 129 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 130 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 131 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 132 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 133 END_2D 134 END DO 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 ) 135 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 ) 136 122 ! … … 145 131 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 146 132 END WHERE 147 END DO 148 DO jl = 1, jpl 149 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 150 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 151 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 152 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 153 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 154 END_3D 155 END DO 156 DO jl = 1, jpl 157 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 158 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 159 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 160 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 161 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 162 END_3D 163 END DO 164 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 165 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 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 ) 166 138 ! 167 139 ! … … 179 151 ENDIF 180 152 zdt = rDt_ice / REAL(icycle) 153 z1_dt = 1._wp / zdt 181 154 182 155 ! --- transport --- ! … … 207 180 !---------------! 208 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 ) 209 188 210 189 ! record at_i before advection (for open water) … … 377 356 ENDIF 378 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 ) 379 371 ! 380 372 !== Open water area ==! … … 384 376 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 385 377 END_2D 386 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1.0_wp ) 387 ! 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 388 388 ! 389 389 ! --- Ensure non-negative fields and in-bound thicknesses --- ! … … 445 445 !! work on H (and not V). It is partly related to the multi-category approach 446 446 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 447 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 448 !! since sv_i and e_i are still good. 447 !! concentration is small). We also limit S and T. 449 448 !!---------------------------------------------------------------------- 450 449 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 490 489 IF( pamsk == 0._wp ) THEN 491 490 DO jl = 1, jpl 492 DO_2D( 1, 0, 1, 0 )491 DO_2D( 0, 0, 1, 0 ) 493 492 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 494 493 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) … … 499 498 ENDIF 500 499 ! 500 END_2D 501 DO_2D( 1, 0, 0, 0 ) 501 502 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 502 503 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) … … 533 534 IF( PRESENT( pua_ho ) ) THEN 534 535 DO jl = 1, jpl 535 DO_2D( 1, 0, 1, 0 ) 536 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 537 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) 538 543 END_2D 539 544 END DO … … 549 554 END_2D 550 555 END DO 551 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1.0_wp )552 556 ! 553 557 END SUBROUTINE adv_umx … … 588 592 ! 589 593 DO jl = 1, jpl !-- flux in x-direction 590 DO_2D( 1, 0, 1, 0 )594 DO_2D( 1, 1, 1, 0 ) 591 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) 592 596 END_2D … … 594 598 ! 595 599 DO jl = 1, jpl !-- first guess of tracer from u-flux 596 DO_2D( 0, 0, 0, 0 )600 DO_2D( 1, 1, 0, 0 ) 597 601 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 598 602 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 601 605 END_2D 602 606 END DO 603 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )604 607 ! 605 608 DO jl = 1, jpl !-- flux in y-direction 606 DO_2D( 1, 0, 1, 0 )609 DO_2D( 1, 0, 0, 0 ) 607 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) 608 611 END_2D … … 612 615 ! 613 616 DO jl = 1, jpl !-- flux in y-direction 614 DO_2D( 1, 0, 1, 0)617 DO_2D( 1, 0, 1, 1 ) 615 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) 616 619 END_2D … … 618 621 ! 619 622 DO jl = 1, jpl !-- first guess of tracer from v-flux 620 DO_2D( 0, 0, 0, 0)623 DO_2D( 0, 0, 1, 1 ) 621 624 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 622 625 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 625 628 END_2D 626 629 END DO 627 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )628 630 ! 629 631 DO jl = 1, jpl !-- flux in x-direction 630 DO_2D( 1, 0, 1, 0 )632 DO_2D( 0, 0, 1, 0 ) 631 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) 632 634 END_2D … … 677 679 ! 678 680 DO jl = 1, jpl 679 DO_2D( 1, 0, 1, 0 )681 DO_2D( 1, 1, 1, 0 ) 680 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 ) 681 685 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 682 686 END_2D … … 695 699 ! 696 700 DO jl = 1, jpl !-- flux in x-direction 697 DO_2D( 1, 0, 1, 0 )701 DO_2D( 1, 1, 1, 0 ) 698 702 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 699 703 END_2D … … 702 706 703 707 DO jl = 1, jpl !-- first guess of tracer from u-flux 704 DO_2D( 0, 0, 0, 0 )708 DO_2D( 1, 1, 0, 0 ) 705 709 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 706 710 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 709 713 END_2D 710 714 END DO 711 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )712 715 713 716 DO jl = 1, jpl !-- flux in y-direction 714 DO_2D( 1, 0, 1, 0 )717 DO_2D( 1, 0, 0, 0 ) 715 718 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 716 719 END_2D … … 721 724 ! 722 725 DO jl = 1, jpl !-- flux in y-direction 723 DO_2D( 1, 0, 1, 0)726 DO_2D( 1, 0, 1, 1 ) 724 727 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 725 728 END_2D … … 728 731 ! 729 732 DO jl = 1, jpl !-- first guess of tracer from v-flux 730 DO_2D( 0, 0, 0, 0)733 DO_2D( 0, 0, 1, 1 ) 731 734 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 732 735 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 735 738 END_2D 736 739 END DO 737 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )738 740 ! 739 741 DO jl = 1, jpl !-- flux in x-direction 740 DO_2D( 1, 0, 1, 0 )742 DO_2D( 0, 0, 1, 0 ) 741 743 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 742 744 END_2D … … 895 897 ! 896 898 DO jl = 1, jpl 897 DO_2D( 1, 0, 1, 0 )899 DO_2D( 0, 0, 1, 0 ) 898 900 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 899 901 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 904 906 ! 905 907 DO jl = 1, jpl 906 DO_2D( 1, 0, 1, 0 )908 DO_2D( 0, 0, 1, 0 ) 907 909 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 908 910 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 914 916 ! 915 917 DO jl = 1, jpl 916 DO_2D( 1, 0, 1, 0 )918 DO_2D( 0, 0, 1, 0 ) 917 919 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 918 920 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 928 930 ! 929 931 DO jl = 1, jpl 930 DO_2D( 1, 0, 1, 0 )932 DO_2D( 0, 0, 1, 0 ) 931 933 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 932 934 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 942 944 ! 943 945 DO jl = 1, jpl 944 DO_2D( 1, 0, 1, 0 )946 DO_2D( 0, 0, 1, 0 ) 945 947 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 946 948 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 963 965 IF( ll_neg ) THEN 964 966 DO jl = 1, jpl 965 DO_2D( 1, 0, 1, 0 )967 DO_2D( 0, 0, 1, 0 ) 966 968 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 967 969 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 973 975 ! !-- High order flux in i-direction --! 974 976 DO jl = 1, jpl 975 DO_2D( 1, 0, 1, 0 )977 DO_2D( 0, 0, 1, 0 ) 976 978 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 977 979 END_2D … … 1031 1033 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1032 1034 DO jl = 1, jpl 1033 DO_2D( 1, 0, 1, 0 )1035 DO_2D( 1, 0, 0, 0 ) 1034 1036 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1035 1037 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1039 1041 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1040 1042 DO jl = 1, jpl 1041 DO_2D( 1, 0, 1, 0 )1043 DO_2D( 1, 0, 0, 0 ) 1042 1044 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1043 1045 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 1048 1050 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1049 1051 DO jl = 1, jpl 1050 DO_2D( 1, 0, 1, 0 )1052 DO_2D( 1, 0, 0, 0 ) 1051 1053 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1052 1054 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1061 1063 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1062 1064 DO jl = 1, jpl 1063 DO_2D( 1, 0, 1, 0 )1065 DO_2D( 1, 0, 0, 0 ) 1064 1066 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1065 1067 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1074 1076 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1075 1077 DO jl = 1, jpl 1076 DO_2D( 1, 0, 1, 0 )1078 DO_2D( 1, 0, 0, 0 ) 1077 1079 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1078 1080 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1095 1097 IF( ll_neg ) THEN 1096 1098 DO jl = 1, jpl 1097 DO_2D( 1, 0, 1, 0 )1099 DO_2D( 1, 0, 0, 0 ) 1098 1100 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1099 1101 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1105 1107 ! !-- High order flux in j-direction --! 1106 1108 DO jl = 1, jpl 1107 DO_2D( 1, 0, 1, 0 )1109 DO_2D( 1, 0, 0, 0 ) 1108 1110 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1109 1111 END_2D … … 1141 1143 ! -------------------------------------------------- 1142 1144 DO jl = 1, jpl 1143 DO_2D( 1, 0, 1, 0 )1145 DO_2D( 0, 0, 1, 0 ) 1144 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 ) 1145 1149 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1146 1150 END_2D … … 1248 1252 ! --------------------------------- 1249 1253 DO jl = 1, jpl 1250 DO_2D( 1, 0, 1, 0 )1254 DO_2D( 0, 0, 1, 0 ) 1251 1255 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1252 1256 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) … … 1259 1263 END_2D 1260 1264 1261 DO_2D( 1, 0, 1, 0 )1265 DO_2D( 1, 0, 0, 0 ) 1262 1266 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1263 1267 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) … … 1616 1620 END SUBROUTINE Hsnow 1617 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 1618 1672 1619 1673 #else -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icedyn_rdgrft.F90
r13495 r13917 349 349 ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 350 350 apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * & 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar ) 352 352 ELSE 353 353 apartf(ji,jl) = 0._wp -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icedyn_rhg_evp.F90
r13497 r13917 129 129 REAl(wp) :: zbetau, zbetav 130 130 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 131 REAL(wp) :: z delta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars131 REAL(wp) :: zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 132 132 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 133 133 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice … … 138 138 REAL(wp) :: zshear, zdum1, zdum2 139 139 ! 140 REAL(wp), DIMENSION(jpi,jpj) :: z p_delt !P/delta at T points140 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 141 141 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 142 142 ! … … 145 145 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 146 146 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 147 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 147 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 148 148 ! 149 149 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 150 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 150 151 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 151 152 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: … … 170 171 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 171 172 !! --- diags 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 173 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 173 175 !! --- SIMIP diags 174 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 368 370 369 371 END_2D 370 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp ) 371 372 DO_2D( 0, 1, 0, 1 ) ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 ! no vector loop 372 373 DO_2D( 0, 0, 0, 0 ) 373 374 374 375 ! shear**2 at T points (doc eq. A16) … … 390 391 391 392 ! delta at T points 392 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 393 394 ! P/delta at T points 395 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 396 393 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 394 395 END_2D 396 CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) 397 398 ! P/delta at T points 399 DO_2D( 1, 1, 1, 1 ) 400 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 401 END_2D 402 403 DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 404 405 ! divergence at T points (duplication to avoid communications) 406 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 407 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 408 & ) * r1_e1e2t(ji,jj) 409 410 ! tension at T points (duplication to avoid communications) 411 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 412 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 413 & ) * r1_e1e2t(ji,jj) 414 397 415 ! alpha for aEVP 398 416 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m … … 411 429 412 430 ! stress at T points (zkt/=0 if landfast) 413 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta *(1._wp - zkt) ) ) * z1_alph1414 zs2(ji,jj) = ( zs2(ji,jj) *zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2431 zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 432 zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 415 433 416 434 END_2D 417 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp )418 435 419 436 ! Save beta at T-points for further computations … … 711 728 & ) * r1_e1e2t(ji,jj) 712 729 zdt2 = zdt * zdt 730 731 zten_i(ji,jj) = zdt 713 732 714 733 ! shear**2 at T points (doc eq. A16) … … 726 745 727 746 ! delta at T points 728 z delta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )729 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta) ) ! 0 if delta=0730 pdelta_i(ji,jj) = z delta + rn_creepl * rswitch747 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 748 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 749 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 731 750 732 751 END_2D 733 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 752 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 753 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 734 754 735 755 ! --- Store the stress tensor for the next time step --- ! 736 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )737 756 pstress1_i (:,:) = zs1 (:,:) 738 757 pstress2_i (:,:) = zs2 (:,:) … … 763 782 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 764 783 765 ! --- stress tensor--- !766 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN767 ! 768 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )784 ! --- Stress tensor invariants (SIMIP diags) --- ! 785 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 786 ! 787 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 769 788 ! 770 DO_2D( 0, 0, 0, 0 ) 771 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 772 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 773 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 774 775 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 776 777 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 778 779 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 780 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 781 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 782 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 783 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 784 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 785 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 786 END_2D 787 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 788 ! 789 CALL iom_put( 'isig1' , zsig1 ) 790 CALL iom_put( 'isig2' , zsig2 ) 791 CALL iom_put( 'isig3' , zsig3 ) 792 ! 793 ! Stress tensor invariants (normal and shear stress N/m) 794 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 795 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 796 797 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 789 DO_2D( 1, 1, 1, 1 ) 790 791 ! Ice stresses 792 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 793 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 794 ! I know, this can be confusing... 795 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 796 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 797 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 798 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 799 800 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 801 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 802 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 803 804 END_2D 805 ! 806 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 807 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 808 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 809 810 DEALLOCATE ( zsig_I, zsig_II ) 811 812 ENDIF 813 814 ! --- Normalized stress tensor principal components --- ! 815 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 816 ! Recommendation 1 : we use ice strength, not replacement pressure 817 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 818 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 819 ! 820 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 821 ! 822 DO_2D( 1, 1, 1, 1 ) 823 824 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 825 ! and **deformations** at current iterates 826 ! following Lemieux & Dupont (2020) 827 zfac = zp_delt(ji,jj) 828 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 829 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 830 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 831 832 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 833 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 834 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 835 836 ! Normalized principal stresses (used to display the ellipse) 837 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 838 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 839 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 840 END_2D 841 ! 842 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 843 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 844 845 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 846 798 847 ENDIF 799 848 … … 931 980 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 932 981 ! close file 933 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid)982 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 934 983 ENDIF 935 984 -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/iceitd.F90
r13472 r13917 627 627 END_2D 628 628 ! 629 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 630 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )631 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )632 !633 DO ji = 1, npti634 jdonor(ji,jl) = jl635 ! how much of a_i you send in cat sup is somewhat arbitrary636 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 637 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 638 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 639 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 640 !! zdaice(ji,jl) = a_i_1d(ji) 641 !! zdvice(ji,jl) = v_i_1d(ji)642 !!clem: these are from UCL and work ok 643 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp644 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1)) * 0.5_wp645 END DO646 !647 IF( npti > 0 ) THEN629 IF( npti > 0 ) THEN 630 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 631 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 632 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 633 ! 634 DO ji = 1, npti 635 jdonor(ji,jl) = jl 636 ! how much of a_i you send in cat sup is somewhat arbitrary 637 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 638 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 639 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 640 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 641 !! zdaice(ji,jl) = a_i_1d(ji) 642 !! zdvice(ji,jl) = v_i_1d(ji) 643 !!clem: these are from UCL and work ok 644 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 645 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 646 END DO 647 ! 648 648 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl=>jl+1 649 649 ! Reset shift parameters … … 666 666 END_2D 667 667 ! 668 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok669 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok670 DO ji = 1, npti671 jdonor(ji,jl) = jl + 1672 zdaice(ji,jl) = a_i_1d(ji)673 zdvice(ji,jl) = v_i_1d(ji)674 END DO675 !676 668 IF( npti > 0 ) THEN 669 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 670 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 671 DO ji = 1, npti 672 jdonor(ji,jl) = jl + 1 673 zdaice(ji,jl) = a_i_1d(ji) 674 zdvice(ji,jl) = v_i_1d(ji) 675 END DO 676 ! 677 677 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl+1=>jl 678 678 ! Reset shift parameters -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icestp.F90
r13472 r13917 55 55 USE icedyn ! sea-ice: dynamics 56 56 USE icethd ! sea-ice: thermodynamics 57 USE icecor ! sea-ice: corrections58 57 USE iceupdate ! sea-ice: sea surface boundary condition update 59 58 USE icedia ! sea-ice: budget diagnostics … … 86 85 PUBLIC ice_init ! called by sbcmod.F90 87 86 87 !! * Substitutions 88 # include "do_loop_substitute.h90" 88 89 !!---------------------------------------------------------------------- 89 90 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 160 161 IF( ln_icedyn .AND. .NOT.lk_c1d ) & 161 162 & CALL ice_dyn( kt, Kmm ) ! -- Ice dynamics 163 ! 164 CALL diag_trends( 1 ) ! record dyn trends 162 165 ! 163 166 ! !== lateral boundary conditions ==! … … 188 191 IF( ln_icethd ) CALL ice_thd( kt ) ! -- Ice thermodynamics 189 192 ! 190 CALL ice_cor( kt , 2 ) ! -- Corrections 191 ! 193 CALL diag_trends( 2 ) ! record thermo trends 192 194 CALL ice_var_glo2eqv ! necessary calls (at least for coupling) 193 195 CALL ice_var_agg( 2 ) ! necessary calls (at least for coupling) … … 196 198 ! 197 199 IF( ln_icediahsb ) CALL ice_dia( kt ) ! -- Diagnostics outputs 200 ! 201 IF( ln_icediachk ) CALL ice_drift_wri( kt ) ! -- Diagnostics outputs for conservation 198 202 ! 199 203 CALL ice_wri( kt ) ! -- Ice outputs … … 208 212 ! --- Ocean time step --- ! 209 213 !-------------------------! 210 IF( ln_icedyn ) CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )! -- update surface ocean stresses214 CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) ) ! -- update surface ocean stresses 211 215 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 212 216 ! … … 257 261 END WHERE 258 262 ! 263 CALL diag_set0 ! set diag of mass, heat and salt fluxes to 0: needed for Agrif child grids 264 ! 259 265 CALL ice_itd_init ! ice thickness distribution initialization 260 266 ! 261 267 CALL ice_thd_init ! set ice thermodynics parameters (clem: important to call it first for melt ponds) 262 268 ! 263 ! ! Initial sea-ice state 264 CALL ice_istate_init 269 CALL ice_sbc_init ! set ice-ocean and ice-atm. coupling parameters 270 ! 271 CALL ice_istate_init ! Initial sea-ice state 265 272 IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 266 273 CALL ice_rst_read( Kbb, Kmm, Kaa ) ! start from a restart file … … 271 278 CALL ice_var_agg(1) 272 279 ! 273 CALL ice_sbc_init ! set ice-ocean and ice-atm. coupling parameters274 !275 280 CALL ice_dyn_init ! set ice dynamics parameters 276 281 ! … … 280 285 ! 281 286 CALL ice_dia_init ! initialization for diags 287 ! 288 CALL ice_drift_init ! initialization for diags of conservation 282 289 ! 283 290 fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction … … 340 347 ENDIF 341 348 ! 342 IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('par_init: online conservation check does not work with BDY')343 !344 349 rDt_ice = REAL(nn_fsbc) * rn_Dt !--- sea-ice timestep and its inverse 345 350 r1_Dt_ice = 1._wp / rDt_ice … … 391 396 !! of the time step 392 397 !!---------------------------------------------------------------------- 393 INTEGER :: ji, jj ! dummy loop index 394 !!---------------------------------------------------------------------- 395 sfx (:,:) = 0._wp ; 396 sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp 397 sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp 398 sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp 399 sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp 400 sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp 401 ! 402 wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp 403 wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp 404 wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp 405 wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp 406 wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp 407 wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp 408 wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 409 wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 410 wfx_snw_sni(:,:) = 0._wp 411 wfx_pnd(:,:) = 0._wp 412 413 hfx_thd(:,:) = 0._wp ; 414 hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp 415 hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp 416 hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp 417 hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp 418 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 419 hfx_err_dif(:,:) = 0._wp 420 wfx_err_sub(:,:) = 0._wp 421 ! 422 diag_heat(:,:) = 0._wp ; diag_sice(:,:) = 0._wp 423 diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp 424 425 ! SIMIP diagnostics 426 qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes 427 t_si (:,:,:) = rt0 ! temp at the ice-snow interface 428 429 tau_icebfr (:,:) = 0._wp ! landfast ice param only (clem: important to keep the init here) 430 cnd_ice (:,:,:) = 0._wp ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 431 qcn_ice (:,:,:) = 0._wp ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 432 qtr_ice_bot(:,:,:) = 0._wp ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 433 qsb_ice_bot(:,:) = 0._wp ! (needed if ln_icethd=F) 434 ! 435 ! for control checks (ln_icediachk) 436 diag_trp_vi(:,:) = 0._wp ; diag_trp_vs(:,:) = 0._wp 437 diag_trp_ei(:,:) = 0._wp ; diag_trp_es(:,:) = 0._wp 438 diag_trp_sv(:,:) = 0._wp 439 398 INTEGER :: ji, jj, jl ! dummy loop index 399 !!---------------------------------------------------------------------- 400 401 DO_2D( 1, 1, 1, 1 ) 402 sfx (ji,jj) = 0._wp ; 403 sfx_bri(ji,jj) = 0._wp ; sfx_lam(ji,jj) = 0._wp 404 sfx_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp 405 sfx_bog(ji,jj) = 0._wp ; sfx_dyn(ji,jj) = 0._wp 406 sfx_bom(ji,jj) = 0._wp ; sfx_sum(ji,jj) = 0._wp 407 sfx_res(ji,jj) = 0._wp ; sfx_sub(ji,jj) = 0._wp 408 ! 409 wfx_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp 410 wfx_sni(ji,jj) = 0._wp ; wfx_opw(ji,jj) = 0._wp 411 wfx_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp 412 wfx_bom(ji,jj) = 0._wp ; wfx_sum(ji,jj) = 0._wp 413 wfx_res(ji,jj) = 0._wp ; wfx_sub(ji,jj) = 0._wp 414 wfx_spr(ji,jj) = 0._wp ; wfx_lam(ji,jj) = 0._wp 415 wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 416 wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 417 wfx_snw_sni(ji,jj) = 0._wp 418 wfx_pnd(ji,jj) = 0._wp 419 420 hfx_thd(ji,jj) = 0._wp ; 421 hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp 422 hfx_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp 423 hfx_bom(ji,jj) = 0._wp ; hfx_sum(ji,jj) = 0._wp 424 hfx_res(ji,jj) = 0._wp ; hfx_sub(ji,jj) = 0._wp 425 hfx_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp 426 hfx_err_dif(ji,jj) = 0._wp 427 wfx_err_sub(ji,jj) = 0._wp 428 ! 429 diag_heat(ji,jj) = 0._wp ; diag_sice(ji,jj) = 0._wp 430 diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp 431 432 tau_icebfr (ji,jj) = 0._wp ! landfast ice param only (clem: important to keep the init here) 433 qsb_ice_bot(ji,jj) = 0._wp ! (needed if ln_icethd=F) 434 435 fhld(ji,jj) = 0._wp ! needed if ln_icethd=F 436 437 ! for control checks (ln_icediachk) 438 diag_trp_vi(ji,jj) = 0._wp ; diag_trp_vs(ji,jj) = 0._wp 439 diag_trp_ei(ji,jj) = 0._wp ; diag_trp_es(ji,jj) = 0._wp 440 diag_trp_sv(ji,jj) = 0._wp 441 ! 442 diag_adv_mass(ji,jj) = 0._wp 443 diag_adv_salt(ji,jj) = 0._wp 444 diag_adv_heat(ji,jj) = 0._wp 445 END_2D 446 447 DO jl = 1, jpl 448 DO_2D( 1, 1, 1, 1 ) 449 ! SIMIP diagnostics 450 t_si (ji,jj,jl) = rt0 ! temp at the ice-snow interface 451 qcn_ice_bot(ji,jj,jl) = 0._wp 452 qcn_ice_top(ji,jj,jl) = 0._wp ! conductive fluxes 453 cnd_ice (ji,jj,jl) = 0._wp ! effective conductivity at the top of ice/snow (ln_cndflx=T) 454 qcn_ice (ji,jj,jl) = 0._wp ! conductive flux (ln_cndflx=T & ln_cndemule=T) 455 qtr_ice_bot(ji,jj,jl) = 0._wp ! part of solar radiation transmitted through the ice needed at least for outputs 456 END_2D 457 ENDDO 458 440 459 END SUBROUTINE diag_set0 460 461 462 SUBROUTINE diag_trends( kn ) 463 !!---------------------------------------------------------------------- 464 !! *** ROUTINE diag_trends *** 465 !! 466 !! ** purpose : diagnostics of the trends. Used for conservation purposes 467 !! and outputs 468 !!---------------------------------------------------------------------- 469 INTEGER, INTENT(in) :: kn ! 1 = after dyn ; 2 = after thermo 470 !!---------------------------------------------------------------------- 471 ! 472 ! --- trends of heat, salt, mass (used for conservation controls) 473 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 474 ! 475 diag_heat(:,:) = diag_heat(:,:) & 476 & - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 477 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 478 diag_sice(:,:) = diag_sice(:,:) & 479 & + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi 480 diag_vice(:,:) = diag_vice(:,:) & 481 & + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi 482 diag_vsnw(:,:) = diag_vsnw(:,:) & 483 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos 484 ! 485 IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend 486 ! 487 ENDIF 488 ! 489 ! --- trends of concentration (used for simip outputs) 490 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 491 ! 492 diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 493 ! 494 IF( kn == 1 ) CALL iom_put( 'afxdyn' , diag_aice ) ! dyn trend 495 IF( kn == 2 ) CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend 496 IF( kn == 2 ) CALL iom_put( 'afxtot' , diag_aice ) ! total trend 497 ! 498 ENDIF 499 ! 500 END SUBROUTINE diag_trends 441 501 442 502 #else -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icetab.F90
r10069 r13917 40 40 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 41 41 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab2d ! input 2D field 42 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT( 42 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) :: tab1d ! output 1D field 43 43 ! 44 44 INTEGER :: jl, jn, jid, jjd … … 61 61 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 62 62 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: tab2d ! input 2D field 63 REAL(wp), DIMENSION(ndim1d) , INTENT( 63 REAL(wp), DIMENSION(ndim1d) , INTENT(inout) :: tab1d ! output 1D field 64 64 ! 65 65 INTEGER :: jn , jid, jjd … … 80 80 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 81 81 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab1d ! input 1D field 82 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( 82 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) :: tab2d ! output 2D field 83 83 ! 84 84 INTEGER :: jl, jn, jid, jjd … … 101 101 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 102 102 REAL(wp), DIMENSION(ndim1d) , INTENT(in ) :: tab1d ! input 1D field 103 REAL(wp), DIMENSION(jpi,jpj), INTENT( 103 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: tab2d ! output 2D field 104 104 ! 105 105 INTEGER :: jn , jid, jjd -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icethd.F90
r13472 r13917 18 18 USE ice ! sea-ice: variables 19 19 !!gm list trop longue ==>>> why not passage en argument d'appel ? 20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot,sprecip, ln_cpl20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 21 21 USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 22 22 & qml_ice, qcn_ice, qtr_ice_top … … 30 30 USE icethd_pnd ! sea-ice: melt ponds 31 31 USE iceitd ! sea-ice: remapping thickness distribution 32 USE icecor ! sea-ice: corrections 32 33 USE icetab ! sea-ice: 1D <==> 2D transformation 33 34 USE icevar ! sea-ice: operations … … 52 53 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F) 53 54 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 54 LOGICAL :: ln_leadhfx ! 55 LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean 55 56 56 57 !! for convergence tests … … 91 92 ! 92 93 INTEGER :: ji, jj, jk, jl ! dummy loop indices 93 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 94 REAL(wp), PARAMETER :: zfric_umin = 0._wp 95 REAL(wp), PARAMETER :: zch = 0.0057_wp 96 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)94 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 95 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 96 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 97 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 97 98 ! 98 99 !!------------------------------------------------------------------- … … 124 125 & ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 125 126 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 127 zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 128 & ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 126 129 END_2D 127 130 ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean … … 130 133 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 131 134 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 135 zvel(ji,jj) = 0._wp 132 136 END_2D 133 137 ENDIF 134 CALL lbc_lnk ( 'icethd', zfric, 'T',1.0_wp )138 CALL lbc_lnk_multi( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 135 139 ! 136 140 !--------------------------------------------------------------------! … … 140 144 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 141 145 ! 142 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget143 ! ! practically no "direct lateral ablation"144 !145 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean146 ! ! temperature and turbulent mixing (McPhee, 1992)147 !148 146 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 149 147 zqld = tmask(ji,jj,1) * rDt_ice * & … … 151 149 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 152 150 153 ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 151 ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 152 ! (mostly<0 but >0 if supercooling) 154 153 zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 155 154 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0 156 157 ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 155 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0 156 157 ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 158 ! (mostly>0 but <0 if supercooling) 158 159 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 159 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 160 161 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 160 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 161 162 162 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 163 163 ! the freezing point, so that we do not have SST < T_freeze 164 ! This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 165 166 !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 167 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 168 169 ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 170 ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 171 IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 172 IF( ln_leadhfx ) THEN ; fhld(ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 173 ELSE ; fhld(ji,jj) = 0._wp 174 ENDIF 164 ! This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 165 ! The following formulation is ok for both normal conditions and supercooling 166 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 167 168 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 169 ! qlead is the energy received from the atm. in the leads. 170 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2) 171 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2) 172 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 173 ! upper bound for fhld: fhld should be equal to zqld 174 ! but we have to make sure that this heat will not make the sst drop below the freezing point 175 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 176 ! The following formulation is ok for both normal conditions and supercooling 177 fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) & ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 178 & - qsb_ice_bot(ji,jj) ) 175 179 qlead(ji,jj) = 0._wp 176 180 ELSE 177 181 fhld (ji,jj) = 0._wp 182 ! upper bound for qlead: qlead should be equal to zqld 183 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 184 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 185 ! and freezing point is reached if zqfr = zqld - qsb*a/dt 186 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 187 ! The following formulation is ok for both normal conditions and supercooling 188 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 178 189 ENDIF 179 190 ! 180 ! Net heat flux on top of the ice-ocean [W.m-2] 181 ! --------------------------------------------- 182 qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 191 ! If ice is landfast and ice concentration reaches its max 192 ! => stop ice formation in open water 193 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp 194 ! 195 ! If the grid cell is almost fully covered by ice (no leads) 196 ! => stop ice formation in open water 197 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp 198 ! 199 ! If ln_leadhfx is false 200 ! => do not use energy of the leads to melt sea-ice 201 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp 202 ! 183 203 END_2D 184 204 … … 191 211 ENDIF 192 212 193 ! ---------------------------------------------------------------------194 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]195 ! ---------------------------------------------------------------------196 ! First step here : non solar + precip - qlead - qsensible197 ! Second step in icethd_dh : heat remaining if total melt (zq_rema)198 ! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar199 qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean200 & - qlead(:,:) * r1_Dt_ice & ! heat flux taken from the ocean where there is open water ice formation201 & - at_i (:,:) * qsb_ice_bot(:,:) & ! heat flux taken by sensible flux202 & - at_i (:,:) * fhld (:,:) ! heat flux taken during bottom growth/melt203 ! ! (fhld should be 0 while bott growth)204 213 !-------------------------------------------------------------------------------------------! 205 214 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 231 240 CALL ice_thd_dh ! Ice-Snow thickness 232 241 CALL ice_thd_pnd ! Melt ponds formation 233 CALL ice_thd_ent( e_i_1d(1:npti,:) , .true.) ! Ice enthalpy remapping242 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 234 243 ENDIF 235 244 CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- ! … … 254 263 ! 255 264 IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- ! 265 ! 266 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 267 ! 268 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice ! ice natural aging incrementation 256 269 ! 257 270 ! convergence tests … … 418 431 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 419 432 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 420 CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )421 433 ! 422 434 ! ocean surface fields 423 435 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 424 436 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 437 CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m ) 425 438 ! 426 439 ! to update ice age … … 510 523 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 511 524 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 512 CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )513 525 ! 514 526 CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) ) -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icethd_dh.F90
r13472 r13917 139 139 ! 140 140 DO ji = 1, npti 141 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 141 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) 142 142 zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 143 143 END DO … … 556 556 ! 557 557 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 558 qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice558 !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 559 559 560 560 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icethd_do.F90
r13295 r13917 131 131 132 132 ! Default new ice thickness 133 WHERE( qlead(:,:) < 0._wp .AND. tau_icebfr(:,:) == 0._wp ) ; ht_i_new(:,:) = rn_hinew ! if cooling and no landfast 134 ELSEWHERE ; ht_i_new(:,:) = 0._wp 133 WHERE( qlead(:,:) < 0._wp ) ! cooling 134 ht_i_new(:,:) = rn_hinew 135 ELSEWHERE 136 ht_i_new(:,:) = 0._wp 135 137 END WHERE 136 138 … … 146 148 ! 147 149 DO_2D( 0, 0, 0, 0 ) 148 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast150 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 149 151 ! -- Wind stress -- ! 150 152 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & … … 198 200 ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 199 201 !------------------------------------------------------------------------------! 200 ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice202 ! it occurs if cooling 201 203 202 204 ! Identify grid points where new ice forms 203 205 npti = 0 ; nptidx(:) = 0 204 206 DO_2D( 1, 1, 1, 1 ) 205 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp) THEN207 IF ( qlead(ji,jj) < 0._wp ) THEN 206 208 npti = npti + 1 207 209 nptidx( npti ) = (jj - 1) * jpi + ji … … 385 387 END DO 386 388 ! --- Ice enthalpy remapping --- ! 387 CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) , .false.)389 CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 388 390 END DO 389 391 -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/icethd_ent.F90
r13472 r13917 38 38 CONTAINS 39 39 40 SUBROUTINE ice_thd_ent( qnew , compute_hfx_err)40 SUBROUTINE ice_thd_ent( qnew ) 41 41 !!------------------------------------------------------------------- 42 42 !! *** ROUTINE ice_thd_ent *** … … 64 64 !!------------------------------------------------------------------- 65 65 REAL(wp), DIMENSION(:,:), INTENT(inout) :: qnew ! new enthlapies (J.m-3, remapped) 66 LOGICAL, INTENT(in) :: compute_hfx_err ! determines whether to compute diag.67 ! error or not68 66 ! 69 67 INTEGER :: ji ! dummy loop indices -
NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/ICE/iceupdate.F90
r13497 r13917 24 24 USE traqsr ! add penetration of solar flux in the calculation of heat budget 25 25 USE icectl ! sea-ice: control prints 26 USE bdy_oce , ONLY : ln_bdy27 26 USE zdfdrg , ONLY : ln_drgice_imp 28 27 ! … … 92 91 ! 93 92 INTEGER :: ji, jj, jl, jk ! dummy loop indices 94 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2)95 93 REAL(wp) :: zqsr ! New solar flux received by the ocean 96 94 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 103 101 WRITE(numout,*)'~~~~~~~~~~~~~~' 104 102 ENDIF 103 104 ! Net heat flux on top of the ice-ocean (W.m-2) 105 !---------------------------------------------- 106 qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 105 107 106 108 ! --- case we bypass ice thermodynamics --- ! … … 115 117 DO_2D( 1, 1, 1, 1 ) 116 118 117 ! Solar heat flux reaching the ocean = zqsr (W.m-2)119 ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2) 118 120 !--------------------------------------------------- 119 121 zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) ) … … 121 123 ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2) 122 124 !--------------------------------------------------- 123 zqmass = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 124 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr 125 126 ! Add the residual from heat diffusion equation and sublimation (W.m-2) 127 !---------------------------------------------------------------------- 128 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) + & 129 & ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 130 125 qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) & 126 & - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) & 127 & + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) & 128 & + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj) 129 131 130 ! New qsr and qns used to compute the oceanic heat flux at the next time step 132 131 !---------------------------------------------------------------------------- 133 qsr(ji,jj) = zqsr 132 ! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice 133 ! else ( cooling or no ice left ), then we suppose that no solar flux has been consumed 134 ! 135 IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN !-- warming and some ice remains 136 ! solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice) 137 qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) & 138 ! + solar flux transmitted thru ice and the 1st ocean level (also not used by sea-ice) 139 & + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) * ( 1._wp - frq_m(ji,jj) ) 140 ! 141 ELSE !-- cooling or no ice left 142 qsr(ji,jj) = zqsr 143 ENDIF 144 ! 145 ! the non-solar is simply derived from the solar flux 134 146 qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr 135 147 136 148 ! Mass flux at the atm. surface 137 149 !----------------------------------- … … 140 152 ! Mass flux at the ocean surface 141 153 !------------------------------------ 142 ! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED) 143 ! ------------------------------------------------------------------------------------- 144 ! The idea of this approach is that the system that we consider is the ICE-OCEAN system 145 ! Thus FW flux = External ( E-P+snow melt) 146 ! Salt flux = Exchanges in the ice-ocean system then converted into FW 147 ! Associated to Ice formation AND Ice melting 148 ! Even if i see Ice melting as a FW and SALT flux 149 ! 150 ! mass flux from ice/ocean 154 ! ice-ocean mass flux 151 155 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 152 156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 153 154 ! add the snow melt water to snow mass flux to the ocean157 158 ! snw-ocean mass flux 155 159 wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) 156 157 ! mass flux at the ocean/ice interface 158 fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) ! F/M mass flux save at least for biogeochemical model 159 emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 160 160 161 ! total mass flux at the ocean/ice interface 162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model 163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux 161 164 162 165 ! Salt flux at the ocean surface … … 262 265 CALL iom_put ('hfxdif' , hfx_dif ) ! heat flux used for ice temperature change 263 266 CALL iom_put ('hfxsnw' , hfx_snw ) ! heat flux used for snow melt 264 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion (included in qt_oce_ai)267 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion 265 268 266 269 ! heat fluxes associated with mass exchange (freeze/melt/precip...) … … 279 282 !--------- 280 283 #if ! defined key_agrif 281 IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final('iceupdate') ! conservation284 IF( ln_icediachk ) CALL ice_cons_final('iceupdate') ! conservation 282 285 #endif 283 IF( ln_icectl 284 IF( sn_cfctl%l_prtctl 285 IF( ln_timing 286 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 287 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('iceupdate') ! prints 288 IF( ln_timing ) CALL timing_stop ('ice_update') ! timing 286 289 ! 287 290 END SUBROUTINE ice_update_flx
Note: See TracChangeset
for help on using the changeset viewer.
