Changeset 10425 for NEMO/trunk/src/ICE
- Timestamp:
- 2018-12-19T22:54:16+01:00 (5 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 19 edited
-
ice.F90 (modified) (2 diffs)
-
ice1d.F90 (modified) (1 diff)
-
icecor.F90 (modified) (1 diff)
-
icectl.F90 (modified) (6 diffs)
-
icedia.F90 (modified) (6 diffs)
-
icedyn.F90 (modified) (3 diffs)
-
icedyn_adv_pra.F90 (modified) (4 diffs)
-
icedyn_adv_umx.F90 (modified) (42 diffs)
-
icedyn_rdgrft.F90 (modified) (4 diffs)
-
icedyn_rhg_evp.F90 (modified) (16 diffs)
-
iceforcing.F90 (modified) (1 diff)
-
iceistate.F90 (modified) (2 diffs)
-
icerst.F90 (modified) (7 diffs)
-
icestp.F90 (modified) (1 diff)
-
icethd.F90 (modified) (3 diffs)
-
icethd_do.F90 (modified) (1 diff)
-
icethd_zdf_bl99.F90 (modified) (13 diffs)
-
iceupdate.F90 (modified) (3 diffs)
-
icewri.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/ice.F90
r10415 r10425 207 207 REAL(wp), PUBLIC, PARAMETER :: epsi20 = 1.e-20_wp !: small number 208 208 209 ! !!** some other parameters for advection using the ULTIMATE-MACHO scheme 210 LOGICAL, PUBLIC, DIMENSION(2) :: l_split_advumx = .FALSE. ! force one iteration at the first time-step 209 211 210 212 ! !!** define arrays … … 463 465 464 466 ice_alloc = MAXVAL( ierr(:) ) 465 IF( ice_alloc /= 0 ) CALL ctl_ warn('ice_alloc: failed to allocate arrays.')467 IF( ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' ) 466 468 ! 467 469 END FUNCTION ice_alloc -
NEMO/trunk/src/ICE/ice1d.F90
r10069 r10425 229 229 230 230 ice1D_alloc = MAXVAL( ierr(:) ) 231 IF( ice1D_alloc /= 0 ) CALL ctl_ warn( 'ice1D_alloc: failed to allocate arrays.')231 IF( ice1D_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice1D_alloc: failed to allocate arrays.' ) 232 232 ! 233 233 END FUNCTION ice1D_alloc -
NEMO/trunk/src/ICE/icecor.F90
r10069 r10425 119 119 END DO 120 120 END DO 121 CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. ) ! lateral boundary conditions121 CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1., v_ice, 'V', -1. ) ! lateral boundary conditions 122 122 ENDIF 123 123 -
NEMO/trunk/src/ICE/icectl.F90
r10415 r10425 77 77 IF( icount == 0 ) THEN 78 78 ! ! water flux 79 pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 79 pdiag_fv = glob_sum( 'icectl', & 80 & -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 80 81 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + & 81 82 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + & … … 84 85 ! 85 86 ! ! salt flux 86 pdiag_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 87 pdiag_fs = glob_sum( 'icectl', & 88 & ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 87 89 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 88 90 & ) * e1e2t(:,:) ) * zconv 89 91 ! 90 92 ! ! heat flux 91 pdiag_ft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 93 pdiag_ft = glob_sum( 'icectl', & 94 & ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 92 95 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 93 96 & ) * e1e2t(:,:) ) * zconv 94 97 95 pdiag_v = glob_sum( SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv )96 97 pdiag_s = glob_sum( SUM( sv_i * rhoi , dim=3 ) * e1e2t * zconv )98 99 pdiag_t = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) &98 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv ) 99 100 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t * zconv ) 101 102 pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 100 103 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 101 104 … … 103 106 104 107 ! water flux 105 zfv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 108 zfv = glob_sum( 'icectl', & 109 & -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 106 110 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + & 107 111 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + & … … 110 114 111 115 ! salt flux 112 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 116 zfs = glob_sum( 'icectl', & 117 & ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 113 118 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 114 119 & ) * e1e2t(:,:) ) * zconv - pdiag_fs 115 120 116 121 ! heat flux 117 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 122 zft = glob_sum( 'icectl', & 123 & ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 118 124 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 119 125 & ) * e1e2t(:,:) ) * zconv - pdiag_ft 120 126 121 127 ! outputs 122 zv = ( ( glob_sum( SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv &128 zv = ( ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv & 123 129 & - pdiag_v ) * r1_rdtice - zfv ) * rday 124 130 125 zs = ( ( glob_sum( SUM( sv_i * rhoi , dim=3 ) * e1e2t ) * zconv &131 zs = ( ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) * zconv & 126 132 & - pdiag_s ) * r1_rdtice + zfs ) * rday 127 133 128 zt = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 134 zt = ( glob_sum( 'icectl', & 135 & ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 129 136 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv & 130 137 & - pdiag_t ) * r1_rdtice + zft 131 138 132 139 ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 133 zvtrp = glob_sum( ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) * zconv * rday134 zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e1e2t ) * zconv135 136 zvmin = glob_min( v_i )137 zamax = glob_max( SUM( a_i, dim=3 ) )138 zamin = glob_min( a_i )140 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) * zconv * rday 141 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) * zconv 142 143 zvmin = glob_min( 'icectl', v_i ) 144 zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 145 zamin = glob_min( 'icectl', a_i ) 139 146 140 147 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 141 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2148 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 142 149 zv_sill = zarea * 2.5e-5 143 150 zs_sill = zarea * 25.e-5 … … 184 191 185 192 ! water flux 186 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday193 zvfx = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 187 194 188 195 ! salt flux 189 zsfx = glob_sum( ( sfx + diag_sice ) * e1e2t ) * zconv * rday196 zsfx = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) * zconv * rday 190 197 191 198 ! heat flux 192 zhfx = glob_sum( ( qt_atm_oi - qt_oce_ai - diag_heat ) * e1e2t ) * zconv199 zhfx = glob_sum( 'icectl', ( qt_atm_oi - qt_oce_ai - diag_heat ) * e1e2t ) * zconv 193 200 194 201 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 195 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2202 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 196 203 zv_sill = zarea * 2.5e-5 197 204 zs_sill = zarea * 25.e-5 … … 400 407 IF( lk_mpp ) THEN 401 408 DO ialert_id = 1, inb_altests 402 CALL mpp_sum( inb_alp(ialert_id))409 CALL mpp_sum('icectl', inb_alp(ialert_id)) 403 410 END DO 404 411 ENDIF -
NEMO/trunk/src/ICE/icedia.F90
r10069 r10425 52 52 ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc ) 53 53 54 IF( lk_mpp ) CALL mpp_sum (ice_dia_alloc )55 IF( ice_dia_alloc /= 0 ) CALL ctl_ warn( 'ice_dia_alloc: failed to allocate arrays')54 CALL mpp_sum ( 'icedia', ice_dia_alloc ) 55 IF( ice_dia_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dia_alloc: failed to allocate arrays' ) 56 56 ! 57 57 END FUNCTION ice_dia_alloc … … 85 85 ! 1 - Contents ! 86 86 ! ----------------------- ! 87 zbg_ivol = glob_sum( vt_i(:,:) * e1e2t(:,:) ) * 1.e-9 ! ice volume (km3)88 zbg_svol = glob_sum( vt_s(:,:) * e1e2t(:,:) ) * 1.e-9 ! snow volume (km3)89 zbg_area = glob_sum( at_i(:,:) * e1e2t(:,:) ) * 1.e-6 ! area (km2)90 zbg_isal = glob_sum( SUM( sv_i(:,:,:), dim=3 ) * e1e2t(:,:) ) * 1.e-9! salt content (pss*km3)91 zbg_item = glob_sum( et_i * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)92 zbg_stem = glob_sum( et_s * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)87 zbg_ivol = glob_sum( 'icedia', vt_i(:,:) * e1e2t(:,:) ) * 1.e-9 ! ice volume (km3) 88 zbg_svol = glob_sum( 'icedia', vt_s(:,:) * e1e2t(:,:) ) * 1.e-9 ! snow volume (km3) 89 zbg_area = glob_sum( 'icedia', at_i(:,:) * e1e2t(:,:) ) * 1.e-6 ! area (km2) 90 zbg_isal = glob_sum( 'icedia', SUM( sv_i(:,:,:), dim=3 ) * e1e2t(:,:) ) * 1.e-9 ! salt content (pss*km3) 91 zbg_item = glob_sum( 'icedia', et_i * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J) 92 zbg_stem = glob_sum( 'icedia', et_s * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J) 93 93 94 94 ! ---------------------------! 95 95 ! 2 - Trends due to forcing ! 96 96 ! ---------------------------! 97 z_frc_volbot = r1_rau0 * glob_sum( -( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-ocean98 z_frc_voltop = r1_rau0 * glob_sum( -( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-atm99 z_frc_sal = r1_rau0 * glob_sum( -sfx(:,:) * e1e2t(:,:) ) * 1.e-9 ! salt fluxes ice/snow-ocean100 z_frc_tembot = glob_sum( qt_oce_ai(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat on top of ocean (and below ice)101 z_frc_temtop = glob_sum( qt_atm_oi(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat on top of ice-coean97 z_frc_volbot = r1_rau0 * glob_sum( 'icedia', -( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-ocean 98 z_frc_voltop = r1_rau0 * glob_sum( 'icedia', -( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-atm 99 z_frc_sal = r1_rau0 * glob_sum( 'icedia', - sfx(:,:) * e1e2t(:,:) ) * 1.e-9 ! salt fluxes ice/snow-ocean 100 z_frc_tembot = glob_sum( 'icedia', qt_oce_ai(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat on top of ocean (and below ice) 101 z_frc_temtop = glob_sum( 'icedia', qt_atm_oi(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat on top of ice-coean 102 102 ! 103 103 frc_voltop = frc_voltop + z_frc_voltop * rdt_ice ! km3 … … 110 110 ! 3 - Content variations ! 111 111 ! ----------------------- ! 112 zdiff_vol = r1_rau0 * glob_sum( ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)113 zdiff_sal = r1_rau0 * glob_sum( ( rhoi* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss)114 zdiff_tem = glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J)112 zdiff_vol = r1_rau0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3) 113 zdiff_sal = r1_rau0 * glob_sum( 'icedia', ( rhoi* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss) 114 zdiff_tem = glob_sum( 'icedia', ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J) 115 115 ! + SUM( qevap_ice * a_i_b, dim=3 ) !! clem: I think this term should not be there (but needs a check) 116 116 … … 125 125 ! 5 - Diagnostics writing ! 126 126 ! ----------------------- ! 127 !!gm I don't understand the division by the ocean surface (i.e. glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt )127 !!gm I don't understand the division by the ocean surface (i.e. glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt ) 128 128 !! and its multiplication bu kt ! is it really what we want ? what is this quantity ? 129 129 !! IF it is really what we want, compute it at kt=nit000, not 3 time by time-step ! … … 135 135 IF( iom_use('ibgheatco') ) CALL iom_put( 'ibgheatco' , zdiff_tem ) ! ice/snow heat content drift (1.e20 J) 136 136 IF( iom_use('ibgheatfx') ) CALL iom_put( 'ibgheatfx' , & ! ice/snow heat flux drift (W/m2) 137 & zdiff_tem /glob_sum( e1e2t(:,:) * 1.e-20 * kt*rdt ) )137 & zdiff_tem /glob_sum( 'icedia', e1e2t(:,:) * 1.e-20 * kt*rdt ) ) 138 138 139 139 IF( iom_use('ibgfrcvoltop') ) CALL iom_put( 'ibgfrcvoltop' , frc_voltop ) ! vol forcing ice/snw-atm (km3 equivalent ocean water) … … 143 143 IF( iom_use('ibgfrctembot') ) CALL iom_put( 'ibgfrctembot' , frc_tembot ) ! heat on top of ocean(below ice) (1.e20 J) 144 144 IF( iom_use('ibgfrchfxtop') ) CALL iom_put( 'ibgfrchfxtop' , & ! heat on top of ice/snw/ocean (W/m2) 145 & frc_temtop / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt )145 & frc_temtop / glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt ) 146 146 IF( iom_use('ibgfrchfxbot') ) CALL iom_put( 'ibgfrchfxbot' , & ! heat on top of ocean(below ice) (W/m2) 147 & frc_tembot / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt )147 & frc_tembot / glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt ) 148 148 149 149 IF( iom_use('ibgvol_tot' ) ) CALL iom_put( 'ibgvol_tot' , zbg_ivol ) ! ice volume (km3) -
NEMO/trunk/src/ICE/icedyn.F90
r10415 r10425 116 116 END DO 117 117 END DO 118 CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )118 CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 119 119 ! 120 120 ! … … 155 155 END DO 156 156 END DO 157 CALL lbc_lnk( zdivu_i, 'T', 1. )157 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 158 158 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 159 159 … … 176 176 END DO 177 177 END DO 178 CALL lbc_lnk( zdivu_i, 'T', 1. )178 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 179 179 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 180 180 -
NEMO/trunk/src/ICE/icedyn_adv_pra.F90
r10069 r10425 101 101 zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 102 102 zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 103 IF( lk_mpp ) CALL mpp_max(zcfl )103 CALL mpp_max( 'icedyn_adv_pra', zcfl ) 104 104 105 105 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp … … 425 425 426 426 !-- Lateral boundary conditions 427 CALL lbc_lnk_multi( psm , 'T', 1., ps0 , 'T', 1. &427 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T', 1., ps0 , 'T', 1. & 428 428 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 429 429 & , psxx, 'T', 1., psyy, 'T', 1. & … … 599 599 600 600 !-- Lateral boundary conditions 601 CALL lbc_lnk_multi( psm , 'T', 1., ps0 , 'T', 1. &601 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T', 1., ps0 , 'T', 1. & 602 602 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 603 603 & , psxx, 'T', 1., psyy, 'T', 1. & … … 640 640 & STAT = ierr ) 641 641 ! 642 IF( lk_mpp ) CALL mpp_sum(ierr )642 CALL mpp_sum( 'icedyn_adv_pra', ierr ) 643 643 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme') 644 644 ! -
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r10418 r10425 35 35 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 36 36 37 REAL(wp), DIMENSION(:,: ), ALLOCATABLE :: amaxu, amaxv37 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: amaxu, amaxv 38 38 39 39 ! advect H all the way (and get V=H*A at the end) … … 119 119 INTEGER :: icycle ! number of sub-timestep for the advection 120 120 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 121 REAL(wp) :: zcfl , zdt 122 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho 123 REAL(wp), DIMENSION(jpi,jpj) :: zhvar 124 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2, z1_ai, z1_aip 121 REAL(wp) :: zdt 122 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv 123 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx 124 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 125 126 127 128 REAL(wp), DIMENSION(jpi,jpj) :: zcu_box, zcv_box 129 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho 130 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip 131 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar 132 125 133 !!---------------------------------------------------------------------- 126 134 ! … … 128 136 ! 129 137 ! 130 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 131 zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 132 zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 133 IF( lk_mpp ) CALL mpp_max( zcfl ) 134 135 IF( zcfl > 0.5 ) THEN ; icycle = 2 136 ELSE ; icycle = 1 138 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 139 ! When needed, the advection split is applied at the next time-step in order to avoid blocking global comm. 140 ! ...this should not affect too much the stability... Was ok on the tests we did... 141 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 142 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 143 144 ! non-blocking global communication send zcflnow and receive zcflprv 145 CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 146 147 IF( zcflprv(1) > .5 ) THEN ; icycle = 2 148 ELSE ; icycle = 1 137 149 ENDIF 138 150 … … 159 171 160 172 IF( ll_zeroup2 ) THEN 161 IF(.NOT. ALLOCATED(amaxu)) ALLOCATE(amaxu (jpi,jpj ))162 IF(.NOT. ALLOCATED(amaxv)) ALLOCATE(amaxv (jpi,jpj ))173 IF(.NOT. ALLOCATED(amaxu)) ALLOCATE(amaxu (jpi,jpj,jpl)) 174 IF(.NOT. ALLOCATED(amaxv)) ALLOCATE(amaxv (jpi,jpj,jpl)) 163 175 ENDIF 164 176 !---------------! … … 167 179 DO jt = 1, icycle 168 180 169 IF( ll_ADVopw ) THEN170 zamsk = 1._wp171 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area172 zamsk = 0._wp173 ELSE181 !!$ IF( ll_ADVopw ) THEN 182 !!$ zamsk = 1._wp 183 !!$ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area 184 !!$ zamsk = 0._wp 185 !!$ ELSE 174 186 zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 175 ENDIF187 !!$ ENDIF 176 188 177 DO jl = 1, jpl 189 WHERE( pa_i(:,:,:) >= epsi20 ) ; z1_ai(:,:,:) = 1._wp / pa_i(:,:,:) 190 ELSEWHERE ; z1_ai(:,:,:) = 0. 191 END WHERE 178 192 ! 179 WHERE( pa_i(:,:,jl) >= epsi20 ) ; z1_ai(:,:) = 1._wp / pa_i(:,:,jl) 180 ELSEWHERE ; z1_ai(:,:) = 0. 181 END WHERE 182 ! 183 WHERE( pa_ip(:,:,jl) >= epsi20 ) ; z1_aip(:,:) = 1._wp / pa_ip(:,:,jl) 184 ELSEWHERE ; z1_aip(:,:) = 0. 185 END WHERE 186 ! 187 IF( ll_zeroup2 ) THEN 193 WHERE( pa_ip(:,:,:) >= epsi20 ) ; z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:) 194 ELSEWHERE ; z1_aip(:,:,:) = 0. 195 END WHERE 196 ! 197 IF( ll_zeroup2 ) THEN 198 DO jl = 1, jpl 188 199 DO jj = 2, jpjm1 189 200 DO ji = fs_2, fs_jpim1 190 amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 191 & pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 192 amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 193 & pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 194 END DO 195 END DO 196 CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.) 197 ENDIF 201 amaxu(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 202 & pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 203 amaxv(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 204 & pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 205 END DO 206 END DO 207 END DO 208 CALL lbc_lnk_multi('icedyn_adv_umx', amaxu, 'T', 1., amaxv, 'T', 1.) 209 ENDIF 210 ! 211 DO jl = 1, jpl 212 zua_ho(:,:,jl) = zudy(:,:) 213 zva_ho(:,:,jl) = zvdx(:,:) 214 END DO 215 216 zamsk = 1._wp 217 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area 218 zamsk = 0._wp 219 ! 220 IF( ll_thickness ) THEN 221 DO jl = 1, jpl 222 zua_ho(:,:,jl) = zudy(:,:) 223 zva_ho(:,:,jl) = zvdx(:,:) 224 END DO 225 ENDIF 198 226 ! 227 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 228 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_i ) ! Ice volume 229 IF( ll_thickness ) pv_i(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 230 ! 231 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 232 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_s ) ! Snw volume 233 IF( ll_thickness ) pv_s(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 234 ! 235 zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 236 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, psv_i ) ! Salt content 237 ! 238 zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 239 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, poa_i ) ! Age content 240 ! 241 DO jk = 1, nlay_i 242 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 243 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content 244 END DO 245 ! 246 DO jk = 1, nlay_s 247 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 248 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content 249 END DO 250 ! 251 IF ( ln_pnd_H12 ) THEN 252 ! 253 DO jl = 1, jpl 254 zua_ho(:,:,jl) = zudy(:,:) 255 zva_ho(:,:,jl) = zvdx(:,:) 256 END DO 257 199 258 zamsk = 1._wp 200 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), & ! Ice area 201 & zua_ho, zva_ho ) 259 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction 202 260 zamsk = 0._wp 203 261 ! 204 IF( ll_thickness ) THEN 205 zua_ho(:,:) = zudy(:,:) 206 zva_ho(:,:) = zvdx(:,:) 207 ENDIF 208 ! 209 zhvar(:,:) = pv_i(:,:,jl) * z1_ai(:,:) 210 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_i (:,:,jl) ) ! Ice volume 211 IF( ll_thickness ) pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 212 ! 213 zhvar(:,:) = pv_s(:,:,jl) * z1_ai(:,:) 214 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_s (:,:,jl) ) ! Snw volume 215 IF( ll_thickness ) pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 216 ! 217 zhvar(:,:) = psv_i(:,:,jl) * z1_ai(:,:) 218 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), psv_i(:,:,jl) ) ! Salt content 219 ! 220 zhvar(:,:) = poa_i(:,:,jl) * z1_ai(:,:) 221 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), poa_i(:,:,jl) ) ! Age content 222 ! 223 DO jk = 1, nlay_i 224 zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai(:,:) 225 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_i(:,:,jk,jl) ) ! Ice heat content 226 END DO 227 ! 228 DO jk = 1, nlay_s 229 zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai(:,:) 230 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_s(:,:,jk,jl) ) ! Snw heat content 231 END DO 232 ! 233 IF ( ln_pnd_H12 ) THEN 234 ! 235 zamsk = 1._wp 236 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), & ! mp fraction 237 & zua_ho, zva_ho ) 238 zamsk = 0._wp 239 240 zhvar(:,:) = pv_ip(:,:,jl) * z1_ai(:,:) 241 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_ip(:,:,jl) ) ! mp volume 242 ENDIF 243 ! 244 ! 245 END DO 246 ! 247 IF( .NOT. ll_ADVopw ) THEN 262 zhvar(:,:,:) = pv_ip(:,:,:) * z1_ai(:,:,:) 263 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume 264 ENDIF 265 ! 266 ! 267 !!$ IF( .NOT. ll_ADVopw ) THEN 248 268 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 249 269 DO jj = 2, jpjm1 … … 253 273 END DO 254 274 END DO 255 CALL lbc_lnk( pato_i(:,:), 'T', 1. )256 ENDIF275 CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. ) 276 !!$ ENDIF 257 277 ! 258 278 END DO … … 279 299 INTEGER , INTENT(in ) :: kt ! number of iteration 280 300 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 281 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2282 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u283 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity284 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: pt ! tracer field285 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: ptc ! tracer content field286 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes287 ! 288 INTEGER :: ji, jj 301 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 302 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 303 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pubox, pvbox ! upstream velocity 304 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pt ! tracer field 305 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptc ! tracer content field 306 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 307 ! 308 INTEGER :: ji, jj, jl ! dummy loop indices 289 309 REAL(wp) :: ztra ! local scalar 290 310 INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 291 REAL(wp), DIMENSION(jpi,jpj ) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt292 REAL(wp), DIMENSION(jpi,jpj ) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc311 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt 312 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc 293 313 !!---------------------------------------------------------------------- 294 314 ! … … 297 317 298 318 IF( ll_gurvan .AND. pamsk==0. ) THEN 299 DO jj = 2, jpjm1 300 DO ji = fs_2, fs_jpim1 301 pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) & 302 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 303 END DO 304 END DO 305 CALL lbc_lnk( pt, 'T', 1. ) 319 DO jl = 1, jpl 320 DO jj = 2, jpjm1 321 DO ji = fs_2, fs_jpim1 322 pt(ji,jj,jl) = ( pt (ji,jj,jl) + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) & 323 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) & 324 & * tmask(ji,jj,1) 325 END DO 326 END DO 327 END DO 328 CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. ) 306 329 ENDIF 307 330 … … 310 333 311 334 ! fluxes in both x-y directions 312 DO jj = 1, jpjm1 313 DO ji = 1, fs_jpim1 314 IF( ll_clem ) THEN 315 zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 316 zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 317 ELSE 318 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 319 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 320 ENDIF 321 END DO 322 END DO 323 324 ELSE 325 ! 326 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 327 ! flux in x-direction 335 DO jl = 1, jpl 328 336 DO jj = 1, jpjm1 329 337 DO ji = 1, fs_jpim1 330 338 IF( ll_clem ) THEN 331 zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 339 zfu_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) 340 zfv_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) 332 341 ELSE 333 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 334 ENDIF 342 zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 343 zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 344 ENDIF 345 END DO 346 END DO 347 END DO 348 349 ELSE 350 ! 351 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 352 ! flux in x-direction 353 DO jl = 1, jpl 354 DO jj = 1, jpjm1 355 DO ji = 1, fs_jpim1 356 IF( ll_clem ) THEN 357 zfu_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) 358 ELSE 359 zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 360 ENDIF 361 END DO 335 362 END DO 336 363 END DO 337 364 338 365 ! first guess of tracer content from u-flux 339 DO jj = 2, jpjm1 340 DO ji = fs_2, fs_jpim1 ! vector opt. 341 IF( ll_clem ) THEN 342 IF( ll_gurvan ) THEN 343 zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 366 DO jl = 1, jpl 367 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 IF( ll_clem ) THEN 370 IF( ll_gurvan ) THEN 371 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 372 ELSE 373 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 374 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 375 & ) * tmask(ji,jj,1) 376 ENDIF 344 377 ELSE 345 zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 346 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 347 & ) * tmask(ji,jj,1) 348 ENDIF 349 ELSE 350 zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) & 351 & * tmask(ji,jj,1) 352 ENDIF 353 !! IF( ji==26 .AND. jj==86) THEN 354 !! WRITE(numout,*) '************************' 355 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 356 !! ENDIF 357 END DO 358 END DO 359 CALL lbc_lnk( zpt, 'T', 1. ) 378 zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 379 & * tmask(ji,jj,1) 380 ENDIF 381 !! IF( ji==26 .AND. jj==86) THEN 382 !! WRITE(numout,*) '************************' 383 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 384 !! ENDIF 385 END DO 386 END DO 387 END DO 388 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 360 389 ! 361 390 ! flux in y-direction 362 DO jj = 1, jpjm1 363 DO ji = 1, fs_jpim1 364 IF( ll_clem ) THEN 365 zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1) 366 ELSE 367 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) 368 ENDIF 369 END DO 370 END DO 371 391 DO jl = 1, jpl 392 DO jj = 1, jpjm1 393 DO ji = 1, fs_jpim1 394 IF( ll_clem ) THEN 395 zfv_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) 396 ELSE 397 zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj+1,jl) 398 ENDIF 399 END DO 400 END DO 401 END DO 372 402 ! 373 403 ELSE !== even ice time step: adv_y then adv_x ==! 374 404 ! flux in y-direction 375 DO jj = 1, jpjm1 376 DO ji = 1, fs_jpim1 377 IF( ll_clem ) THEN 378 zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 379 ELSE 380 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 381 ENDIF 405 DO jl = 1, jpl 406 DO jj = 1, jpjm1 407 DO ji = 1, fs_jpim1 408 IF( ll_clem ) THEN 409 zfv_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) 410 ELSE 411 zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 412 ENDIF 413 END DO 382 414 END DO 383 415 END DO 384 416 385 417 ! first guess of tracer content from v-flux 386 DO jj = 2, jpjm1 387 DO ji = fs_2, fs_jpim1 ! vector opt. 388 IF( ll_clem ) THEN 389 IF( ll_gurvan ) THEN 390 zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 418 DO jl = 1, jpl 419 DO jj = 2, jpjm1 420 DO ji = fs_2, fs_jpim1 ! vector opt. 421 IF( ll_clem ) THEN 422 IF( ll_gurvan ) THEN 423 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 424 ELSE 425 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 426 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 427 & * tmask(ji,jj,1) 428 ENDIF 391 429 ELSE 392 zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 393 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 394 & * tmask(ji,jj,1) 395 ENDIF 396 ELSE 397 zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 398 & * tmask(ji,jj,1) 399 ENDIF 400 !! IF( ji==26 .AND. jj==86) THEN 401 !! WRITE(numout,*) '************************' 402 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 403 !! ENDIF 404 END DO 405 END DO 406 CALL lbc_lnk( zpt, 'T', 1. ) 430 zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 431 & * tmask(ji,jj,1) 432 ENDIF 433 !! IF( ji==26 .AND. jj==86) THEN 434 !! WRITE(numout,*) '************************' 435 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 436 !! ENDIF 437 END DO 438 END DO 439 END DO 440 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 407 441 ! 408 442 ! flux in x-direction 409 DO jj = 1, jpjm1 410 DO ji = 1, fs_jpim1 411 IF( ll_clem ) THEN 412 zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj) 413 ELSE 414 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) 415 ENDIF 443 DO jl = 1, jpl 444 DO jj = 1, jpjm1 445 DO ji = 1, fs_jpim1 446 IF( ll_clem ) THEN 447 zfu_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) 448 ELSE 449 zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * zpt(ji+1,jj,jl) 450 ENDIF 451 END DO 416 452 END DO 417 453 END DO … … 425 461 426 462 IF( ll_zeroup2 ) THEN 427 DO jj = 1, jpjm1 428 DO ji = 1, fs_jpim1 ! vector opt. 429 IF( amaxu(ji,jj) == 0._wp ) zfu_ups(ji,jj) = 0._wp 430 IF( amaxv(ji,jj) == 0._wp ) zfv_ups(ji,jj) = 0._wp 463 DO jl = 1, jpl 464 DO jj = 1, jpjm1 465 DO ji = 1, fs_jpim1 ! vector opt. 466 IF( amaxu(ji,jj,jl) == 0._wp ) zfu_ups(ji,jj,jl) = 0._wp 467 IF( amaxv(ji,jj,jl) == 0._wp ) zfv_ups(ji,jj,jl) = 0._wp 468 END DO 431 469 END DO 432 470 END DO … … 434 472 435 473 ! guess after content field with upstream scheme 436 DO jj = 2, jpjm1 437 DO ji = fs_2, fs_jpim1 438 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj ) & 439 & + zfv_ups(ji,jj) - zfv_ups(ji ,jj-1) ) * r1_e1e2t(ji,jj) 440 IF( ll_clem ) THEN 441 IF( ll_gurvan ) THEN 442 zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 474 DO jl = 1, jpl 475 DO jj = 2, jpjm1 476 DO ji = fs_2, fs_jpim1 477 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj ,jl) & 478 & + zfv_ups(ji,jj,jl) - zfv_ups(ji ,jj-1,jl) ) * r1_e1e2t(ji,jj) 479 IF( ll_clem ) THEN 480 IF( ll_gurvan ) THEN 481 zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 482 ELSE 483 zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) & 484 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 485 & * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 486 ENDIF 443 487 ELSE 444 zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra + ( pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) & 445 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 446 & * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 488 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 447 489 ENDIF 448 ELSE 449 zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 450 ENDIF 451 !! IF( ji==26 .AND. jj==86) THEN 452 !! WRITE(numout,*) '**************************' 453 !! WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 454 !! ENDIF 490 !! IF( ji==26 .AND. jj==86) THEN 491 !! WRITE(numout,*) '**************************' 492 !! WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 493 !! ENDIF 494 END DO 455 495 END DO 456 496 END DO 457 CALL lbc_lnk( zt_ups, 'T', 1. )497 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 458 498 459 499 ! High order (_ho) fluxes … … 476 516 ! final trend with corrected fluxes 477 517 ! ------------------------------------ 478 DO jj = 2, jpjm1 479 DO ji = fs_2, fs_jpim1 480 IF( ll_gurvan ) THEN 481 ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 482 ELSE 483 ztra = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) & 484 & + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 485 & + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 486 ENDIF 487 pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 488 489 IF( pt(ji,jj) < -epsi20 ) THEN 490 WRITE(numout,*) 'T<0 ',pt(ji,jj) 491 ENDIF 492 493 IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 ) pt(ji,jj) = 0._wp 494 495 !! IF( ji==26 .AND. jj==86) THEN 496 !! WRITE(numout,*) 'zt high order',pt(ji,jj) 497 !! ENDIF 498 END DO 499 END DO 500 CALL lbc_lnk( pt, 'T', 1. ) 518 DO jl = 1, jpl 519 DO jj = 2, jpjm1 520 DO ji = fs_2, fs_jpim1 521 IF( ll_gurvan ) THEN 522 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 523 ELSE 524 ztra = ( - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) & 525 & + ( pt(ji,jj,jl) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 526 & + ( pt(ji,jj,jl) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 527 ENDIF 528 pt(ji,jj,jl) = ( pt(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 529 530 IF( pt(ji,jj,jl) < -epsi20 ) THEN 531 WRITE(numout,*) 'T<0 ',pt(ji,jj,jl) 532 ENDIF 533 534 IF( pt(ji,jj,jl) < 0._wp .AND. pt(ji,jj,jl) >= -epsi20 ) pt(ji,jj,jl) = 0._wp 535 536 !! IF( ji==26 .AND. jj==86) THEN 537 !! WRITE(numout,*) 'zt high order',pt(ji,jj) 538 !! ENDIF 539 END DO 540 END DO 541 END DO 542 CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. ) 501 543 ENDIF 502 544 … … 505 547 IF( ll_clem ) THEN 506 548 IF( pamsk == 0. ) THEN 507 DO jj = 1, jpjm1 508 DO ji = 1, fs_jpim1 509 IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 510 zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj) 511 zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj) 512 ELSE 513 zfu_ho (ji,jj) = 0._wp 514 zfu_ups(ji,jj) = 0._wp 515 ENDIF 516 ! 517 IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 518 zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj) 519 zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj) 520 ELSE 521 zfv_ho (ji,jj) = 0._wp 522 zfv_ups(ji,jj) = 0._wp 523 ENDIF 524 ENDDO 525 ENDDO 549 DO jl = 1, jpl 550 DO jj = 1, jpjm1 551 DO ji = 1, fs_jpim1 552 IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 553 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 554 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 555 ELSE 556 zfu_ho (ji,jj,jl) = 0._wp 557 zfu_ups(ji,jj,jl) = 0._wp 558 ENDIF 559 ! 560 IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 561 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 562 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 563 ELSE 564 zfv_ho (ji,jj,jl) = 0._wp 565 zfv_ups(ji,jj,jl) = 0._wp 566 ENDIF 567 END DO 568 END DO 569 END DO 526 570 ENDIF 527 571 ENDIF 528 572 529 573 IF( ll_zeroup5 ) THEN 530 DO jj = 2, jpjm1 531 DO ji = 2, fs_jpim1 ! vector opt. 532 zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 533 & - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 534 IF( zpt(ji,jj) < 0. ) THEN 535 zfu_ho(ji,jj) = zfu_ups(ji,jj) 536 zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj) 537 zfv_ho(ji,jj) = zfv_ups(ji,jj) 538 zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1) 539 ENDIF 540 END DO 541 END DO 542 CALL lbc_lnk_multi( zfu_ho, 'U', -1., zfv_ho, 'V', -1. ) 574 DO jl = 1, jpl 575 DO jj = 2, jpjm1 576 DO ji = 2, fs_jpim1 ! vector opt. 577 zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 578 & - ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 579 IF( zpt(ji,jj,jl) < 0. ) THEN 580 zfu_ho(ji ,jj,jl) = zfu_ups(ji ,jj,jl) 581 zfu_ho(ji-1,jj,jl) = zfu_ups(ji-1,jj,jl) 582 zfv_ho(ji ,jj,jl) = zfv_ups(ji ,jj,jl) 583 zfv_ho(ji,jj-1,jl) = zfv_ups(ji,jj-1,jl) 584 ENDIF 585 END DO 586 END DO 587 END DO 588 CALL lbc_lnk_multi( 'icedyn_adv_umx', zfu_ho, 'U', -1., zfv_ho, 'V', -1. ) 543 589 ENDIF 544 590 … … 546 592 ! ---------------------------- 547 593 IF( PRESENT( pua_ho ) ) THEN 548 DO jj = 1, jpjm1 549 DO ji = 1, fs_jpim1 550 pua_ho(ji,jj) = zfu_ho(ji,jj) 551 pva_ho(ji,jj) = zfv_ho(ji,jj) 594 DO jl = 1, jpl 595 DO jj = 1, jpjm1 596 DO ji = 1, fs_jpim1 597 pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) 598 pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) 599 END DO 552 600 END DO 553 601 END DO … … 558 606 ! final trend with corrected fluxes 559 607 ! ------------------------------------ 560 DO jj = 2, jpjm1 561 DO ji = fs_2, fs_jpim1 562 ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) * pdt 563 564 ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1) 565 566 !! IF( ji==26 .AND. jj==86) THEN 567 !! WRITE(numout,*) 'ztc high order',ptc(ji,jj) 568 !! ENDIF 569 570 END DO 571 END DO 572 CALL lbc_lnk( ptc, 'T', 1. ) 608 DO jl = 1, jpl 609 DO jj = 2, jpjm1 610 DO ji = fs_2, fs_jpim1 611 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt 612 613 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1) 614 615 !! IF( ji==26 .AND. jj==86) THEN 616 !! WRITE(numout,*) 'ztc high order',ptc(ji,jj) 617 !! ENDIF 618 619 END DO 620 END DO 621 END DO 622 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. ) 573 623 ENDIF 574 624 … … 593 643 INTEGER , INTENT(in ) :: kt ! number of iteration 594 644 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 595 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pt ! tracer fields596 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components597 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components598 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step599 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes600 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content601 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes602 ! 603 INTEGER :: ji, jj ! dummy loop indices645 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields 646 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu, pv ! 2 ice velocity components 647 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 648 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptc ! tracer content at before time step 649 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 650 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt_ups ! upstream guess of tracer content 651 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes 652 ! 653 INTEGER :: ji, jj, jl ! dummy loop indices 604 654 LOGICAL :: ll_xy = .TRUE. 605 REAL(wp), DIMENSION(jpi,jpj ) :: zzt655 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zzt 606 656 !!---------------------------------------------------------------------- 607 657 ! 608 658 IF( .NOT.ll_xy ) THEN !-- no alternate directions --! 609 659 ! 610 DO jj = 1, jpjm1 611 DO ji = 1, fs_jpim1 612 IF( ll_clem ) THEN 613 pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 614 pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 615 ELSE 616 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 617 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 618 ENDIF 660 DO jl = 1, jpl 661 DO jj = 1, jpjm1 662 DO ji = 1, fs_jpim1 663 IF( ll_clem ) THEN 664 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 665 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 666 ELSE 667 pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 668 pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 669 ENDIF 670 END DO 619 671 END DO 620 672 END DO … … 638 690 ! 639 691 ! flux in x-direction 640 DO jj = 1, jpjm1 641 DO ji = 1, fs_jpim1 642 IF( ll_clem ) THEN 643 pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 644 ELSE 645 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 646 ENDIF 692 DO jl = 1, jpl 693 DO jj = 1, jpjm1 694 DO ji = 1, fs_jpim1 695 IF( ll_clem ) THEN 696 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 697 ELSE 698 pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 699 ENDIF 700 END DO 647 701 END DO 648 702 END DO … … 651 705 652 706 ! first guess of tracer content from u-flux 653 DO jj = 2, jpjm1 654 DO ji = fs_2, fs_jpim1 ! vector opt. 655 IF( ll_clem ) THEN 656 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 657 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 707 DO jl = 1, jpl 708 DO jj = 2, jpjm1 709 DO ji = fs_2, fs_jpim1 ! vector opt. 710 IF( ll_clem ) THEN 711 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 712 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 658 713 & * tmask(ji,jj,1) 659 ELSE 660 zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 661 ENDIF 662 END DO 663 END DO 664 CALL lbc_lnk( zzt, 'T', 1. ) 714 ELSE 715 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 716 ENDIF 717 END DO 718 END DO 719 END DO 720 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 665 721 666 722 ! flux in y-direction 667 DO jj = 1, jpjm1 668 DO ji = 1, fs_jpim1 669 IF( ll_clem ) THEN 670 pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) 671 ELSE 672 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) 673 ENDIF 723 DO jl = 1, jpl 724 DO jj = 1, jpjm1 725 DO ji = 1, fs_jpim1 726 IF( ll_clem ) THEN 727 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 728 ELSE 729 pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 730 ENDIF 731 END DO 674 732 END DO 675 733 END DO … … 680 738 ! 681 739 ! flux in y-direction 682 DO jj = 1, jpjm1 683 DO ji = 1, fs_jpim1 684 IF( ll_clem ) THEN 685 pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 686 ELSE 687 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 688 ENDIF 740 DO jl = 1, jpl 741 DO jj = 1, jpjm1 742 DO ji = 1, fs_jpim1 743 IF( ll_clem ) THEN 744 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 745 ELSE 746 pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 747 ENDIF 748 END DO 689 749 END DO 690 750 END DO … … 693 753 ! 694 754 ! first guess of tracer content from v-flux 695 DO jj = 2, jpjm1 696 DO ji = fs_2, fs_jpim1 ! vector opt. 697 IF( ll_clem ) THEN 698 zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 699 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 700 & * tmask(ji,jj,1) 701 ELSE 702 zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 703 ENDIF 704 END DO 705 END DO 706 CALL lbc_lnk( zzt, 'T', 1. ) 755 DO jl = 1, jpl 756 DO jj = 2, jpjm1 757 DO ji = fs_2, fs_jpim1 ! vector opt. 758 IF( ll_clem ) THEN 759 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 760 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 761 & * tmask(ji,jj,1) 762 ELSE 763 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 764 ENDIF 765 END DO 766 END DO 767 END DO 768 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 707 769 ! 708 770 ! flux in x-direction 709 DO jj = 1, jpjm1 710 DO ji = 1, fs_jpim1 711 IF( ll_clem ) THEN 712 pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) 713 ELSE 714 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) 715 ENDIF 771 DO jl = 1, jpl 772 DO jj = 1, jpjm1 773 DO ji = 1, fs_jpim1 774 IF( ll_clem ) THEN 775 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 776 ELSE 777 pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 778 ENDIF 779 END DO 716 780 END DO 717 781 END DO … … 749 813 INTEGER , INTENT(in ) :: kt ! number of iteration 750 814 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 751 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pt ! tracer fields752 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components753 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components754 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity755 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step756 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points757 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes758 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content759 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes760 ! 761 INTEGER :: ji, jj ! dummy loop indices815 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields 816 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu, pv ! 2 ice velocity components 817 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 818 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pubox, pvbox ! upstream velocity 819 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptc ! tracer content at before time step 820 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 821 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 822 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt_ups ! upstream guess of tracer content 823 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes 824 ! 825 INTEGER :: ji, jj, jl ! dummy loop indices 762 826 REAL(wp) :: ztra 763 REAL(wp), DIMENSION(jpi,jpj ) :: zzt, zzfu_ho, zzfv_ho827 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zzt, zzfu_ho, zzfv_ho 764 828 !!---------------------------------------------------------------------- 765 829 ! … … 776 840 777 841 ! first guess of tracer content from u-flux 778 DO jj = 2, jpjm1 779 DO ji = fs_2, fs_jpim1 ! vector opt. 780 IF( ll_clem ) THEN 842 DO jl = 1, jpl 843 DO jj = 2, jpjm1 844 DO ji = fs_2, fs_jpim1 ! vector opt. 845 IF( ll_clem ) THEN 846 IF( ll_gurvan ) THEN 847 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 848 ELSE 849 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 850 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 851 & ) * tmask(ji,jj,1) 852 ENDIF 853 ELSE 854 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 855 ENDIF 856 END DO 857 END DO 858 END DO 859 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 860 861 ELSE 862 863 DO jl = 1, jpl 864 DO jj = 2, jpjm1 865 DO ji = fs_2, fs_jpim1 ! vector opt. 781 866 IF( ll_gurvan ) THEN 782 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 867 zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) & 868 & - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) 783 869 ELSE 784 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 785 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 786 & ) * tmask(ji,jj,1) 787 ENDIF 788 ELSE 789 zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 790 ENDIF 791 END DO 792 END DO 793 CALL lbc_lnk( zzt, 'T', 1. ) 794 795 ELSE 796 797 DO jj = 2, jpjm1 798 DO ji = fs_2, fs_jpim1 ! vector opt. 799 IF( ll_gurvan ) THEN 800 zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 801 & - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) 802 ELSE 803 zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 804 & - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 805 ENDIF 806 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 807 END DO 808 END DO 809 CALL lbc_lnk( zzt, 'T', 1. ) 870 zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) & 871 & - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 872 ENDIF 873 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 874 END DO 875 END DO 876 END DO 877 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 810 878 ENDIF 811 879 ! … … 832 900 833 901 ! first guess of tracer content from v-flux 834 DO jj = 2, jpjm1 835 DO ji = fs_2, fs_jpim1 ! vector opt. 836 IF( ll_clem ) THEN 837 IF( ll_gurvan ) THEN 838 zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 902 DO jl = 1, jpl 903 DO jj = 2, jpjm1 904 DO ji = fs_2, fs_jpim1 ! vector opt. 905 IF( ll_clem ) THEN 906 IF( ll_gurvan ) THEN 907 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 908 ELSE 909 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 910 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 911 & ) * tmask(ji,jj,1) 912 ENDIF 839 913 ELSE 840 zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 841 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 842 & ) * tmask(ji,jj,1) 843 ENDIF 844 ELSE 845 zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 846 & * tmask(ji,jj,1) 847 ENDIF 848 END DO 849 END DO 850 CALL lbc_lnk( zzt, 'T', 1. ) 914 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 915 & * tmask(ji,jj,1) 916 ENDIF 917 END DO 918 END DO 919 END DO 920 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 851 921 852 922 ELSE 853 923 854 DO jj = 2, jpjm1 855 DO ji = fs_2, fs_jpim1 856 IF( ll_gurvan ) THEN 857 zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 858 & - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) 859 ELSE 860 zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 861 & - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 862 ENDIF 863 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 864 END DO 865 END DO 866 CALL lbc_lnk( zzt, 'T', 1. ) 924 DO jl = 1, jpl 925 DO jj = 2, jpjm1 926 DO ji = fs_2, fs_jpim1 927 IF( ll_gurvan ) THEN 928 zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) & 929 & - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) 930 ELSE 931 zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) & 932 & - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 933 ENDIF 934 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 935 END DO 936 END DO 937 END DO 938 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 867 939 ENDIF 868 940 ! … … 889 961 ENDIF 890 962 ELSE 891 zzfu_ho(:,: ) = pfu_ho(:,:)892 zzfv_ho(:,: ) = pfv_ho(:,:)963 zzfu_ho(:,:,:) = pfu_ho(:,:,:) 964 zzfv_ho(:,:,:) = pfv_ho(:,:,:) 893 965 ! 1st iteration of nonosc (limit the flux with the upstream solution) 894 966 IF( ll_clem ) THEN … … 898 970 ENDIF 899 971 ! guess after content field with high order 900 DO jj = 2, jpjm1 901 DO ji = fs_2, fs_jpim1 902 ztra = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 903 zzt(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 904 END DO 905 END DO 906 CALL lbc_lnk( zzt, 'T', 1. ) 972 DO jl = 1, jpl 973 DO jj = 2, jpjm1 974 DO ji = fs_2, fs_jpim1 975 ztra = - ( zzfu_ho(ji,jj,jl) - zzfu_ho(ji-1,jj,jl) + zzfv_ho(ji,jj,jl) - zzfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 976 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 977 END DO 978 END DO 979 END DO 980 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 907 981 ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 908 982 IF( ll_clem ) THEN … … 930 1004 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 931 1005 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 932 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pu ! ice i-velocity component933 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: puc ! ice i-velocity * A component934 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pt ! tracer fields935 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: pt_u ! tracer at u-point936 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: pfu_ho ! high order flux937 ! 938 INTEGER :: ji, jj ! dummy loop indices1006 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu ! ice i-velocity component 1007 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc ! ice i-velocity * A component 1008 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields 1009 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u ! tracer at u-point 1010 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho ! high order flux 1011 ! 1012 INTEGER :: ji, jj, jl ! dummy loop indices 939 1013 REAL(wp) :: zcu, zdx2, zdx4 ! - - 940 REAL(wp), DIMENSION(jpi,jpj ) :: ztu1, ztu2, ztu3, ztu41014 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4 941 1015 !!---------------------------------------------------------------------- 942 1016 ! 943 1017 ! !-- Laplacian in i-direction --! 944 DO jj = 2, jpjm1 ! First derivative (gradient) 945 DO ji = 1, fs_jpim1 946 ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 947 END DO 948 ! ! Second derivative (Laplacian) 949 DO ji = fs_2, fs_jpim1 950 ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 1018 DO jl = 1, jpl 1019 DO jj = 2, jpjm1 ! First derivative (gradient) 1020 DO ji = 1, fs_jpim1 1021 ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 1022 END DO 1023 ! ! Second derivative (Laplacian) 1024 DO ji = fs_2, fs_jpim1 1025 ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 1026 END DO 951 1027 END DO 952 1028 END DO 953 CALL lbc_lnk( ztu2, 'T', 1. )1029 CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. ) 954 1030 ! 955 1031 ! !-- BiLaplacian in i-direction --! 956 DO jj = 2, jpjm1 ! Third derivative 957 DO ji = 1, fs_jpim1 958 ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 959 END DO 960 ! ! Fourth derivative 961 DO ji = fs_2, fs_jpim1 962 ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 1032 DO jl = 1, jpl 1033 DO jj = 2, jpjm1 ! Third derivative 1034 DO ji = 1, fs_jpim1 1035 ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 1036 END DO 1037 ! ! Fourth derivative 1038 DO ji = fs_2, fs_jpim1 1039 ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 1040 END DO 963 1041 END DO 964 1042 END DO 965 CALL lbc_lnk( ztu4, 'T', 1. )1043 CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. ) 966 1044 ! 967 1045 ! … … 970 1048 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 971 1049 ! 972 DO jj = 1, jpjm1 973 DO ji = 1, fs_jpim1 ! vector opt. 974 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 975 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 1050 DO jl = 1, jpl 1051 DO jj = 1, jpjm1 1052 DO ji = 1, fs_jpim1 ! vector opt. 1053 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 1054 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 1055 END DO 976 1056 END DO 977 1057 END DO … … 979 1059 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 980 1060 ! 981 DO jj = 1, jpjm1 982 DO ji = 1, fs_jpim1 ! vector opt. 983 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 984 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 985 & - zcu * ( pt(ji+1,jj) - pt(ji,jj) ) ) 1061 DO jl = 1, jpl 1062 DO jj = 1, jpjm1 1063 DO ji = 1, fs_jpim1 ! vector opt. 1064 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1065 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 1066 & - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 1067 END DO 986 1068 END DO 987 1069 END DO … … 989 1071 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 990 1072 ! 991 DO jj = 1, jpjm1 992 DO ji = 1, fs_jpim1 ! vector opt. 993 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 994 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1073 DO jl = 1, jpl 1074 DO jj = 1, jpjm1 1075 DO ji = 1, fs_jpim1 ! vector opt. 1076 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1077 zdx2 = e1u(ji,jj) * e1u(ji,jj) 995 1078 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 996 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 997 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 998 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 999 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) 1079 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 1080 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 1081 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 1082 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 1083 END DO 1000 1084 END DO 1001 1085 END DO … … 1003 1087 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1004 1088 ! 1005 DO jj = 1, jpjm1 1006 DO ji = 1, fs_jpim1 ! vector opt. 1007 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1008 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1089 DO jl = 1, jpl 1090 DO jj = 1, jpjm1 1091 DO ji = 1, fs_jpim1 ! vector opt. 1092 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1093 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1009 1094 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 1010 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 1011 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 1012 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 1013 & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) 1095 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 1096 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 1097 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 1098 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 1099 END DO 1014 1100 END DO 1015 1101 END DO … … 1017 1103 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1018 1104 ! 1019 DO jj = 1, jpjm1 1020 DO ji = 1, fs_jpim1 ! vector opt. 1021 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1022 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1023 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 1024 zdx4 = zdx2 * zdx2 1025 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 1026 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 1027 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 1028 & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) & 1029 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj) & 1030 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) ) 1105 DO jl = 1, jpl 1106 DO jj = 1, jpjm1 1107 DO ji = 1, fs_jpim1 ! vector opt. 1108 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1109 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1110 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 1111 zdx4 = zdx2 * zdx2 1112 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 1113 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 1114 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 1115 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 1116 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) & 1117 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 1118 END DO 1031 1119 END DO 1032 1120 END DO … … 1035 1123 ! !-- High order flux in i-direction --! 1036 1124 IF( ll_neg ) THEN 1125 DO jl = 1, jpl 1126 DO jj = 1, jpjm1 1127 DO ji = 1, fs_jpim1 1128 IF( pt_u(ji,jj,jl) < 0._wp ) THEN 1129 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 1130 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 1131 ENDIF 1132 END DO 1133 END DO 1134 END DO 1135 ENDIF 1136 1137 DO jl = 1, jpl 1037 1138 DO jj = 1, jpjm1 1038 DO ji = 1, fs_jpim1 1039 IF( pt_u(ji,jj) < 0._wp ) THEN 1040 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 1041 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 1139 DO ji = 1, fs_jpim1 ! vector opt. 1140 IF( ll_clem ) THEN 1141 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 1142 ELSE 1143 pfu_ho(ji,jj,jl) = puc(ji,jj,jl) * pt_u(ji,jj,jl) 1042 1144 ENDIF 1043 1145 END DO 1044 END DO1045 ENDIF1046 1047 DO jj = 1, jpjm11048 DO ji = 1, fs_jpim1 ! vector opt.1049 IF( ll_clem ) THEN1050 pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj)1051 ELSE1052 pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)1053 ENDIF1054 1146 END DO 1055 1147 END DO … … 1071 1163 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 1072 1164 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1073 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pv ! ice j-velocity component1074 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity*A component1075 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pt ! tracer fields1076 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: pt_v ! tracer at v-point1077 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: pfv_ho ! high order flux1078 ! 1079 INTEGER :: ji, jj ! dummy loop indices1165 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pv ! ice j-velocity component 1166 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pvc ! ice j-velocity*A component 1167 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields 1168 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_v ! tracer at v-point 1169 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfv_ho ! high order flux 1170 ! 1171 INTEGER :: ji, jj, jl ! dummy loop indices 1080 1172 REAL(wp) :: zcv, zdy2, zdy4 ! - - 1081 REAL(wp), DIMENSION(jpi,jpj ) :: ztv1, ztv2, ztv3, ztv41173 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4 1082 1174 !!---------------------------------------------------------------------- 1083 1175 ! 1084 1176 ! !-- Laplacian in j-direction --! 1085 DO jj = 1, jpjm1 ! First derivative (gradient) 1086 DO ji = fs_2, fs_jpim1 1087 ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1177 DO jl = 1, jpl 1178 DO jj = 1, jpjm1 ! First derivative (gradient) 1179 DO ji = fs_2, fs_jpim1 1180 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1181 END DO 1182 END DO 1183 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 1184 DO ji = fs_2, fs_jpim1 1185 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1186 END DO 1088 1187 END DO 1089 1188 END DO 1090 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 1091 DO ji = fs_2, fs_jpim1 1092 ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 1189 CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. ) 1190 ! 1191 ! !-- BiLaplacian in j-direction --! 1192 DO jl = 1, jpl 1193 DO jj = 1, jpjm1 ! First derivative 1194 DO ji = fs_2, fs_jpim1 1195 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1196 END DO 1197 END DO 1198 DO jj = 2, jpjm1 ! Second derivative 1199 DO ji = fs_2, fs_jpim1 1200 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1201 END DO 1093 1202 END DO 1094 1203 END DO 1095 CALL lbc_lnk( ztv2, 'T', 1. ) 1096 ! 1097 ! !-- BiLaplacian in j-direction --! 1098 DO jj = 1, jpjm1 ! First derivative 1099 DO ji = fs_2, fs_jpim1 1100 ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1101 END DO 1102 END DO 1103 DO jj = 2, jpjm1 ! Second derivative 1104 DO ji = fs_2, fs_jpim1 1105 ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 1106 END DO 1107 END DO 1108 CALL lbc_lnk( ztv4, 'T', 1. ) 1204 CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. ) 1109 1205 ! 1110 1206 ! 1111 1207 SELECT CASE (kn_umx ) 1112 !1208 ! 1113 1209 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1114 DO jj = 1, jpjm1 1115 DO ji = 1, fs_jpim1 1116 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 1117 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 1210 DO jl = 1, jpl 1211 DO jj = 1, jpjm1 1212 DO ji = 1, fs_jpim1 1213 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1214 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1215 END DO 1118 1216 END DO 1119 1217 END DO 1120 1218 ! 1121 1219 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1122 DO jj = 1, jpjm1 1123 DO ji = 1, fs_jpim1 1124 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1125 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 1126 & - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) 1127 END DO 1128 END DO 1129 CALL lbc_lnk( pt_v, 'V', 1. ) 1220 DO jl = 1, jpl 1221 DO jj = 1, jpjm1 1222 DO ji = 1, fs_jpim1 1223 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1224 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1225 & - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1226 END DO 1227 END DO 1228 END DO 1229 CALL lbc_lnk( 'icedyn_adv_umx', pt_v, 'V', 1. ) 1130 1230 ! 1131 1231 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1132 DO jj = 1, jpjm1 1133 DO ji = 1, fs_jpim1 1134 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1135 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1232 DO jl = 1, jpl 1233 DO jj = 1, jpjm1 1234 DO ji = 1, fs_jpim1 1235 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1236 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1136 1237 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1137 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 1138 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 1139 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 1140 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 1238 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1239 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1240 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1241 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1242 END DO 1141 1243 END DO 1142 1244 END DO 1143 1245 ! 1144 1246 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1145 DO jj = 1, jpjm1 1146 DO ji = 1, fs_jpim1 1147 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1148 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1247 DO jl = 1, jpl 1248 DO jj = 1, jpjm1 1249 DO ji = 1, fs_jpim1 1250 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1251 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1149 1252 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1150 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 1151 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 1152 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 1153 & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 1253 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1254 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1255 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1256 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1257 END DO 1154 1258 END DO 1155 1259 END DO 1156 1260 ! 1157 1261 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1158 DO jj = 1, jpjm1 1159 DO ji = 1, fs_jpim1 1160 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1161 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1262 DO jl = 1, jpl 1263 DO jj = 1, jpjm1 1264 DO ji = 1, fs_jpim1 1265 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1266 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1162 1267 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1163 zdy4 = zdy2 * zdy2 1164 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 1165 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 1166 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 1167 & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) & 1168 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj) & 1169 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) ) 1268 zdy4 = zdy2 * zdy2 1269 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1270 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1271 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1272 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 1273 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) & 1274 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 1275 END DO 1170 1276 END DO 1171 1277 END DO … … 1174 1280 ! !-- High order flux in j-direction --! 1175 1281 IF( ll_neg ) THEN 1282 DO jl = 1, jpl 1283 DO jj = 1, jpjm1 1284 DO ji = 1, fs_jpim1 1285 IF( pt_v(ji,jj,jl) < 0._wp ) THEN 1286 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1287 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1288 ENDIF 1289 END DO 1290 END DO 1291 END DO 1292 ENDIF 1293 1294 DO jl = 1, jpl 1176 1295 DO jj = 1, jpjm1 1177 DO ji = 1, fs_jpim1 1178 IF( pt_v(ji,jj) < 0._wp ) THEN 1179 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 1180 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 1296 DO ji = 1, fs_jpim1 ! vector opt. 1297 IF( ll_clem ) THEN 1298 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1299 ELSE 1300 pfv_ho(ji,jj,jl) = pvc(ji,jj,jl) * pt_v(ji,jj,jl) 1181 1301 ENDIF 1182 1302 END DO 1183 END DO1184 ENDIF1185 1186 DO jj = 1, jpjm11187 DO ji = 1, fs_jpim1 ! vector opt.1188 IF( ll_clem ) THEN1189 pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj)1190 ELSE1191 pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)1192 ENDIF1193 1303 END DO 1194 1304 END DO … … 1212 1322 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 1213 1323 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1214 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e21215 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a1216 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pv ! ice j-velocity => v*e11217 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity *A => v*e1*a1218 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: ptc, pt, pt_low ! before field & upstream guess of after field1219 REAL(wp), DIMENSION ( jpi,jpj), INTENT(inout) :: pfv_low, pfu_low ! upstream flux1220 REAL(wp), DIMENSION ( jpi,jpj), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux1221 ! 1222 INTEGER :: ji, jj ! dummy loop indices1324 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2 1325 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a 1326 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice j-velocity => v*e1 1327 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pvc ! ice j-velocity *A => v*e1*a 1328 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: ptc, pt, pt_low ! before field & upstream guess of after field 1329 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_low, pfu_low ! upstream flux 1330 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux 1331 ! 1332 INTEGER :: ji, jj, jl ! dummy loop indices 1223 1333 REAL(wp) :: zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2 ! local scalars 1224 1334 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef ! - - 1225 REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt 1335 REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo 1336 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_low, ztj_low, zzt 1226 1337 !!---------------------------------------------------------------------- 1227 1338 zbig = 1.e+40_wp 1228 1339 zsml = epsi20 1229 1340 1230 1341 IF( ll_zeroup2 ) THEN 1231 DO jj = 1, jpjm1 1232 DO ji = 1, fs_jpim1 ! vector opt. 1233 IF( amaxu(ji,jj) == 0._wp ) pfu_ho(ji,jj) = 0._wp 1234 IF( amaxv(ji,jj) == 0._wp ) pfv_ho(ji,jj) = 0._wp 1342 DO jl = 1, jpl 1343 DO jj = 1, jpjm1 1344 DO ji = 1, fs_jpim1 ! vector opt. 1345 IF( amaxu(ji,jj,jl) == 0._wp ) pfu_ho(ji,jj,jl) = 0._wp 1346 IF( amaxv(ji,jj,jl) == 0._wp ) pfv_ho(ji,jj,jl) = 0._wp 1347 END DO 1235 1348 END DO 1236 1349 END DO 1237 1350 ENDIF 1238 1351 1239 1352 IF( ll_zeroup4 ) THEN 1240 DO jj = 1, jpjm1 1241 DO ji = 1, fs_jpim1 ! vector opt. 1242 IF( pfu_low(ji,jj) == 0._wp ) pfu_ho(ji,jj) = 0._wp 1243 IF( pfv_low(ji,jj) == 0._wp ) pfv_ho(ji,jj) = 0._wp 1353 DO jl = 1, jpl 1354 DO jj = 1, jpjm1 1355 DO ji = 1, fs_jpim1 ! vector opt. 1356 IF( pfu_low(ji,jj,jl) == 0._wp ) pfu_ho(ji,jj,jl) = 0._wp 1357 IF( pfv_low(ji,jj,jl) == 0._wp ) pfv_ho(ji,jj,jl) = 0._wp 1358 END DO 1244 1359 END DO 1245 1360 END DO … … 1248 1363 1249 1364 IF( ll_zeroup1 ) THEN 1250 DO jj = 2, jpjm1 1251 DO ji = fs_2, fs_jpim1 1252 IF( ll_gurvan ) THEN 1253 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1254 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1255 ELSE 1256 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1257 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 1258 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1259 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1260 & ) * tmask(ji,jj,1) 1261 ENDIF 1262 IF( zzt(ji,jj) < 0._wp ) THEN 1263 pfu_ho(ji,jj) = pfu_low(ji,jj) 1264 pfv_ho(ji,jj) = pfv_low(ji,jj) 1265 WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj) 1266 ENDIF 1267 !! IF( ji==26 .AND. jj==86) THEN 1268 !! WRITE(numout,*) 'zzt high order',zzt(ji,jj) 1269 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1270 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1271 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 1272 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 1273 !! ENDIF 1274 IF( ll_gurvan ) THEN 1275 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1276 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1277 ELSE 1278 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1279 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 1280 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1281 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1282 & ) * tmask(ji,jj,1) 1283 ENDIF 1284 IF( zzt(ji,jj) < 0._wp ) THEN 1285 pfu_ho(ji-1,jj) = pfu_low(ji-1,jj) 1286 pfv_ho(ji,jj-1) = pfv_low(ji,jj-1) 1287 WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj) 1288 ENDIF 1289 IF( ll_gurvan ) THEN 1290 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1291 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1292 ELSE 1293 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1294 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 1295 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1296 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1297 & ) * tmask(ji,jj,1) 1298 ENDIF 1299 IF( zzt(ji,jj) < 0._wp ) THEN 1300 WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj) 1301 ENDIF 1302 END DO 1303 END DO 1304 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) 1365 DO jl = 1, jpl 1366 DO jj = 2, jpjm1 1367 DO ji = fs_2, fs_jpim1 1368 IF( ll_gurvan ) THEN 1369 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1370 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1371 ELSE 1372 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1373 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 1374 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1375 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1376 & ) * tmask(ji,jj,1) 1377 ENDIF 1378 IF( zzt(ji,jj,jl) < 0._wp ) THEN 1379 pfu_ho(ji,jj,jl) = pfu_low(ji,jj,jl) 1380 pfv_ho(ji,jj,jl) = pfv_low(ji,jj,jl) 1381 WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 1382 ENDIF 1383 !! IF( ji==26 .AND. jj==86) THEN 1384 !! WRITE(numout,*) 'zzt high order',zzt(ji,jj) 1385 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1386 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1387 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1388 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt 1389 !! ENDIF 1390 IF( ll_gurvan ) THEN 1391 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1392 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1393 ELSE 1394 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1395 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 1396 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1397 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1398 & ) * tmask(ji,jj,1) 1399 ENDIF 1400 IF( zzt(ji,jj,jl) < 0._wp ) THEN 1401 pfu_ho(ji-1,jj,jl) = pfu_low(ji-1,jj,jl) 1402 pfv_ho(ji,jj-1,jl) = pfv_low(ji,jj-1,jl) 1403 WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 1404 ENDIF 1405 IF( ll_gurvan ) THEN 1406 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1407 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1408 ELSE 1409 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1410 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 1411 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1412 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1413 & ) * tmask(ji,jj,1) 1414 ENDIF 1415 IF( zzt(ji,jj,jl) < 0._wp ) THEN 1416 WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 1417 ENDIF 1418 END DO 1419 END DO 1420 END DO 1421 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) 1305 1422 ENDIF 1306 1423 … … 1308 1425 ! antidiffusive flux : high order minus low order 1309 1426 ! -------------------------------------------------- 1310 DO jj = 1, jpjm1 1311 DO ji = 1, fs_jpim1 ! vector opt. 1312 pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj) 1313 pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj) 1314 END DO 1427 DO jl = 1, jpl 1428 DO jj = 1, jpjm1 1429 DO ji = 1, fs_jpim1 ! vector opt. 1430 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_low(ji,jj,jl) 1431 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_low(ji,jj,jl) 1432 END DO 1433 END DO 1315 1434 END DO 1316 1435 … … 1325 1444 IF( ll_prelimiter_zalesak ) THEN 1326 1445 1327 DO jj = 2, jpjm1 1328 DO ji = fs_2, fs_jpim1 1329 zti_low(ji,jj)= pt_low(ji+1,jj ) 1330 ztj_low(ji,jj)= pt_low(ji ,jj+1) 1331 END DO 1332 END DO 1333 CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 1446 DO jl = 1, jpl 1447 DO jj = 2, jpjm1 1448 DO ji = fs_2, fs_jpim1 1449 zti_low(ji,jj,jl)= pt_low(ji+1,jj ,jl) 1450 ztj_low(ji,jj,jl)= pt_low(ji ,jj+1,jl) 1451 END DO 1452 END DO 1453 END DO 1454 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. ) 1334 1455 1335 1456 !! this does not work ?? … … 1353 1474 !! END DO 1354 1475 1355 DO jj = 2, jpjm1 1356 DO ji = fs_2, fs_jpim1 1357 IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND. & 1358 & pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN 1359 ! 1360 IF( pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND. & 1361 & pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 1362 ! 1363 IF( pfu_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji-1,jj) ) <= 0 .AND. & 1364 & pfv_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji,jj-1) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 1365 ! 1366 ENDIF 1367 END DO 1368 END DO 1369 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1476 DO jl = 1, jpl 1477 DO jj = 2, jpjm1 1478 DO ji = fs_2, fs_jpim1 1479 IF ( pfu_ho(ji,jj,jl) * ( pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) <= 0. .AND. & 1480 & pfv_ho(ji,jj,jl) * ( pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) <= 0. ) THEN 1481 ! 1482 IF( pfu_ho(ji,jj,jl) * ( zti_low(ji+1,jj,jl) - zti_low(ji,jj,jl) ) <= 0 .AND. & 1483 & pfv_ho(ji,jj,jl) * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj,jl) ) <= 0) pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 1484 ! 1485 IF( pfu_ho(ji,jj,jl) * ( pt_low(ji ,jj,jl) - pt_low(ji-1,jj,jl) ) <= 0 .AND. & 1486 & pfv_ho(ji,jj,jl) * ( pt_low(ji ,jj,jl) - pt_low(ji,jj-1,jl) ) <= 0) pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 1487 ! 1488 ENDIF 1489 END DO 1490 END DO 1491 END DO 1492 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1370 1493 1371 1494 ELSEIF( ll_prelimiter_devore ) THEN 1372 DO jj = 2, jpjm1 1373 DO ji = fs_2, fs_jpim1 1374 zti_low(ji,jj)= pt_low(ji+1,jj ) 1375 ztj_low(ji,jj)= pt_low(ji ,jj+1) 1376 END DO 1377 END DO 1378 CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 1495 DO jl = 1, jpl 1496 DO jj = 2, jpjm1 1497 DO ji = fs_2, fs_jpim1 1498 zti_low(ji,jj,jl)= pt_low(ji+1,jj ,jl) 1499 ztj_low(ji,jj,jl)= pt_low(ji ,jj+1,jl) 1500 END DO 1501 END DO 1502 END DO 1503 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. ) 1379 1504 1380 1505 z1_dt = 1._wp / pdt 1381 DO jj = 2, jpjm1 1382 DO ji = fs_2, fs_jpim1 1383 zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) ) 1384 pfu_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , & 1385 & zsign * ( pt_low (ji ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji ,jj) * z1_dt , & 1386 & zsign * ( zti_low(ji+1,jj) - zti_low(ji ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) ) 1387 1388 zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) ) 1389 pfv_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , & 1390 & zsign * ( pt_low (ji,jj ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj ) * z1_dt , & 1391 & zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj ) ) * e1e2t(ji,jj+1) * z1_dt ) ) 1392 END DO 1393 END DO 1394 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1395 1506 DO jl = 1, jpl 1507 DO jj = 2, jpjm1 1508 DO ji = fs_2, fs_jpim1 1509 zsign = SIGN( 1., pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) 1510 pfu_ho(ji,jj,jl) = zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj,jl)) , & 1511 & zsign * ( pt_low (ji ,jj,jl) - pt_low (ji-1,jj,jl) ) * e1e2t(ji ,jj) * z1_dt , & 1512 & zsign * ( zti_low(ji+1,jj,jl) - zti_low(ji ,jj,jl) ) * e1e2t(ji+1,jj) * z1_dt ) ) 1513 1514 zsign = SIGN( 1., pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) 1515 pfv_ho(ji,jj,jl) = zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj,jl)) , & 1516 & zsign * ( pt_low (ji,jj ,jl) - pt_low (ji,jj-1,jl) ) * e1e2t(ji,jj ) * z1_dt , & 1517 & zsign * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj ,jl) ) * e1e2t(ji,jj+1) * z1_dt ) ) 1518 END DO 1519 END DO 1520 END DO 1521 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1522 1396 1523 ENDIF 1397 1398 1524 1525 1399 1526 ! Search local extrema 1400 1527 ! -------------------- 1401 1528 ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover 1402 DO jj = 1, jpj1403 DO ji = 1, jpi1404 IF ( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN1405 zbup(ji,jj) = -zbig1406 zbdo(ji,jj) = zbig1407 ELSEIF( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) > 0._wp ) THEN1408 zbup(ji,jj) = pt_low(ji,jj)1409 zbdo(ji,jj) = pt_low(ji,jj)1410 ELSEIF( pt(ji,jj) > 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN1411 zbup(ji,jj) = pt(ji,jj)1412 zbdo(ji,jj) = pt(ji,jj)1413 ELSE1414 zbup(ji,jj) = MAX( pt(ji,jj) , pt_low(ji,jj) )1415 zbdo(ji,jj) = MIN( pt(ji,jj) , pt_low(ji,jj) )1416 ENDIF1417 END DO1418 END DO1419 1420 1421 1529 z1_dt = 1._wp / pdt 1422 DO jj = 2, jpjm1 1423 DO ji = fs_2, fs_jpim1 ! vector opt. 1424 ! 1425 IF( .NOT. ll_9points ) THEN 1426 zup = MAX( zbup(ji,jj), zbup(ji-1,jj ), zbup(ji+1,jj ), zbup(ji ,jj-1), zbup(ji ,jj+1) ) ! search max/min in neighbourhood 1427 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 1428 ! 1429 ELSE 1430 zup = MAX( zbup(ji,jj), zbup(ji-1,jj ), zbup(ji+1,jj ), zbup(ji ,jj-1), zbup(ji ,jj+1), & ! search max/min in neighbourhood 1431 & zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1) ) 1432 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1), & 1433 & zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1) ) 1434 ENDIF 1435 ! 1436 zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji ,jj) ) & ! positive/negative part of the flux 1437 & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj ) ) 1438 zneg = MAX( 0., pfu_ho(ji ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & 1439 & + MAX( 0., pfv_ho(ji,jj ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 1440 ! 1441 IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 1442 zneg2 = ( pt(ji,jj) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1443 & ) * ( 1. - pamsk ) 1444 zpos2 = ( - pt(ji,jj) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1445 & ) * ( 1. - pamsk ) 1446 ELSE 1447 zneg2 = 0. ; zpos2 = 0. 1448 ENDIF 1449 ! 1450 ! ! up & down beta terms 1451 IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 1452 ELSE ; zbetup(ji,jj) = 0. ! zbig 1453 ENDIF 1454 ! 1455 IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 1456 ELSE ; zbetdo(ji,jj) = 0. ! zbig 1457 ENDIF 1458 ! 1459 ! if all the points are outside ice cover 1460 IF( zup == -zbig ) zbetup(ji,jj) = 0. ! zbig 1461 IF( zdo == zbig ) zbetdo(ji,jj) = 0. ! zbig 1462 ! 1530 DO jl = 1, jpl 1531 1532 DO jj = 1, jpj 1533 DO ji = 1, jpi 1534 IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN 1535 zbup(ji,jj) = -zbig 1536 zbdo(ji,jj) = zbig 1537 ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) > 0._wp ) THEN 1538 zbup(ji,jj) = pt_low(ji,jj,jl) 1539 zbdo(ji,jj) = pt_low(ji,jj,jl) 1540 ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN 1541 zbup(ji,jj) = pt(ji,jj,jl) 1542 zbdo(ji,jj) = pt(ji,jj,jl) 1543 ELSE 1544 zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_low(ji,jj,jl) ) 1545 zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_low(ji,jj,jl) ) 1546 ENDIF 1547 END DO 1548 END DO 1549 1550 DO jj = 2, jpjm1 1551 DO ji = fs_2, fs_jpim1 ! vector opt. 1552 ! 1553 IF( .NOT. ll_9points ) THEN 1554 zup = MAX( zbup(ji,jj), zbup(ji-1,jj ), zbup(ji+1,jj ), zbup(ji ,jj-1), zbup(ji ,jj+1) ) ! search max/min in neighbourhood 1555 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 1556 ! 1557 ELSE 1558 zup = MAX( zbup(ji,jj), zbup(ji-1,jj ), zbup(ji+1,jj ), zbup(ji ,jj-1), zbup(ji ,jj+1), & ! search max/min in neighbourhood 1559 & zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1) ) 1560 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1), & 1561 & zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1) ) 1562 ENDIF 1563 ! 1564 zpos = MAX( 0., pfu_ho(ji-1,jj,jl) ) - MIN( 0., pfu_ho(ji ,jj,jl) ) & ! positive/negative part of the flux 1565 & + MAX( 0., pfv_ho(ji,jj-1,jl) ) - MIN( 0., pfv_ho(ji,jj ,jl) ) 1566 zneg = MAX( 0., pfu_ho(ji ,jj,jl) ) - MIN( 0., pfu_ho(ji-1,jj,jl) ) & 1567 & + MAX( 0., pfv_ho(ji,jj ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) ) 1568 ! 1569 IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 1570 zneg2 = ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1571 & ) * ( 1. - pamsk ) 1572 zpos2 = ( - pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1573 & ) * ( 1. - pamsk ) 1574 ELSE 1575 zneg2 = 0. ; zpos2 = 0. 1576 ENDIF 1577 ! 1578 ! ! up & down beta terms 1579 IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_low(ji,jj,jl) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 1580 ELSE ; zbetup(ji,jj,jl) = 0. ! zbig 1581 ENDIF 1582 ! 1583 IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_low(ji,jj,jl) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 1584 ELSE ; zbetdo(ji,jj,jl) = 0. ! zbig 1585 ENDIF 1586 ! 1587 ! if all the points are outside ice cover 1588 IF( zup == -zbig ) zbetup(ji,jj,jl) = 0. ! zbig 1589 IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0. ! zbig 1590 ! 1463 1591 1464 1592 !! IF( ji==26 .AND. jj==86) THEN … … 1489 1617 ! ENDIF 1490 1618 ! 1619 END DO 1491 1620 END DO 1492 1621 END DO 1493 CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)1622 CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 1494 1623 1495 1624 1496 1625 ! monotonic flux in the y direction 1497 1626 ! --------------------------------- 1498 DO jj = 1, jpjm1 1499 DO ji = 1, fs_jpim1 ! vector opt. 1500 zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) 1501 zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) 1502 zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj) ) 1503 ! 1504 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1505 1506 pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj) 1627 DO jl = 1, jpl 1628 DO jj = 1, jpjm1 1629 DO ji = 1, fs_jpim1 ! vector opt. 1630 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1631 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 1632 zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj,jl) ) 1633 ! 1634 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1635 1636 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_low(ji,jj,jl) 1507 1637 1508 1638 !! IF( ji==26 .AND. jj==86) THEN 1509 1639 !! WRITE(numout,*) 'coefU',zcoef 1510 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj )) * r1_e1e2t(ji,jj) * pdt1511 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj )) * r1_e1e2t(ji,jj) * pdt1640 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1641 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1512 1642 !! ENDIF 1513 1643 1514 END DO1515 END DO1516 1517 DO jj = 1, jpjm11518 DO ji = 1, fs_jpim1 ! vector opt.1519 zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )1520 zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )1521 zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj) )1522 !1523 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )1524 1525 pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj)1644 END DO 1645 END DO 1646 1647 DO jj = 1, jpjm1 1648 DO ji = 1, fs_jpim1 ! vector opt. 1649 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1650 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 1651 zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj,jl) ) 1652 ! 1653 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 1654 1655 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_low(ji,jj,jl) 1526 1656 1527 1657 !! IF( ji==26 .AND. jj==86) THEN 1528 1658 !! WRITE(numout,*) 'coefV',zcoef 1529 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj )) * r1_e1e2t(ji,jj) * pdt1530 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1 )) * r1_e1e2t(ji,jj) * pdt1659 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1660 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt 1531 1661 !! ENDIF 1532 END DO 1662 END DO 1663 END DO 1664 1665 ! clem test 1666 DO jj = 2, jpjm1 1667 DO ji = 2, fs_jpim1 ! vector opt. 1668 IF( ll_gurvan ) THEN 1669 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1670 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1671 ELSE 1672 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1673 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 1674 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1675 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1676 & ) * tmask(ji,jj,1) 1677 ENDIF 1678 IF( zzt(ji,jj,jl) < -epsi20 ) THEN 1679 WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj,jl) 1680 ENDIF 1681 END DO 1682 END DO 1683 1533 1684 END DO 1534 1685 1535 ! clem test1536 DO jj = 2, jpjm11537 DO ji = 2, fs_jpim1 ! vector opt.1538 IF( ll_gurvan ) THEN1539 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &1540 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)1541 ELSE1542 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &1543 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &1544 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1545 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1546 & ) * tmask(ji,jj,1)1547 ENDIF1548 IF( zzt(ji,jj) < -epsi20 ) THEN1549 WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj)1550 ENDIF1551 END DO1552 END DO1553 1554 1686 ! 1555 1687 ! … … 1562 1694 !! ** Purpose : compute flux limiter 1563 1695 !!---------------------------------------------------------------------- 1564 REAL(wp) 1565 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e21566 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a1567 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pt ! ice tracer1568 REAL(wp), DIMENSION ( jpi,jpj), INTENT(inout) :: pfu_ho ! high order flux1569 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ), OPTIONAL :: pfu_ups ! upstream flux1696 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1697 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2 1698 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a 1699 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer 1700 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfu_ho ! high order flux 1701 REAL(wp), DIMENSION (:,:,:), INTENT(in ), OPTIONAL :: pfu_ups ! upstream flux 1570 1702 ! 1571 1703 REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr 1572 INTEGER :: ji, jj ! dummy loop indices 1573 REAL(wp), DIMENSION (jpi,jpj) :: zslpx ! tracer slopes 1574 !!---------------------------------------------------------------------- 1575 ! 1576 DO jj = 2, jpjm1 1577 DO ji = fs_2, fs_jpim1 ! vector opt. 1578 zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1) 1704 INTEGER :: ji, jj, jl ! dummy loop indices 1705 REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpx ! tracer slopes 1706 !!---------------------------------------------------------------------- 1707 ! 1708 DO jl = 1, jpl 1709 DO jj = 2, jpjm1 1710 DO ji = fs_2, fs_jpim1 ! vector opt. 1711 zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 1712 END DO 1579 1713 END DO 1580 1714 END DO 1581 CALL lbc_lnk( zslpx, 'U', -1.) ! lateral boundary cond.1715 CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond. 1582 1716 1583 DO jj = 2, jpjm1 1584 DO ji = fs_2, fs_jpim1 ! vector opt. 1585 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1586 1587 Rjm = zslpx(ji-1,jj) 1588 Rj = zslpx(ji ,jj) 1589 Rjp = zslpx(ji+1,jj) 1590 1591 IF( PRESENT(pfu_ups) ) THEN 1592 1593 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 1594 ELSE ; Rr = Rjp 1717 DO jl = 1, jpl 1718 DO jj = 2, jpjm1 1719 DO ji = fs_2, fs_jpim1 ! vector opt. 1720 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1721 1722 Rjm = zslpx(ji-1,jj,jl) 1723 Rj = zslpx(ji ,jj,jl) 1724 Rjp = zslpx(ji+1,jj,jl) 1725 1726 IF( PRESENT(pfu_ups) ) THEN 1727 1728 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 1729 ELSE ; Rr = Rjp 1730 ENDIF 1731 1732 zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1733 IF( Rj > 0. ) THEN 1734 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), & 1735 & MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1736 ELSE 1737 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), & 1738 & MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1739 ENDIF 1740 pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 1741 1742 ELSE 1743 IF( Rj /= 0. ) THEN 1744 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1745 ELSE ; Cr = Rjp / Rj 1746 ENDIF 1747 ELSE 1748 Cr = 0. 1749 !IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 1750 !ELSE ; Cr = Rjp * 1.e20 1751 !ENDIF 1752 ENDIF 1753 1754 ! -- superbee -- 1755 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1756 ! -- van albada 2 -- 1757 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1758 1759 ! -- sweby (with beta=1) -- 1760 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1761 ! -- van Leer -- 1762 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1763 ! -- ospre -- 1764 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1765 ! -- koren -- 1766 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1767 ! -- charm -- 1768 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1769 !ELSE ; zpsi = 0. 1770 !ENDIF 1771 ! -- van albada 1 -- 1772 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1773 ! -- smart -- 1774 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1775 ! -- umist -- 1776 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1777 1778 ! high order flux corrected by the limiter 1779 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 1780 1595 1781 ENDIF 1596 1597 zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj) 1598 IF( Rj > 0. ) THEN 1599 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)), & 1600 & MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 1601 ELSE 1602 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)), & 1603 & MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 1604 ENDIF 1605 pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter 1606 1607 ELSE 1608 IF( Rj /= 0. ) THEN 1609 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1610 ELSE ; Cr = Rjp / Rj 1611 ENDIF 1612 ELSE 1613 Cr = 0. 1614 !IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 1615 !ELSE ; Cr = Rjp * 1.e20 1616 !ENDIF 1617 ENDIF 1618 1619 ! -- superbee -- 1620 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1621 ! -- van albada 2 -- 1622 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1623 1624 ! -- sweby (with beta=1) -- 1625 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1626 ! -- van Leer -- 1627 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1628 ! -- ospre -- 1629 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1630 ! -- koren -- 1631 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1632 ! -- charm -- 1633 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1634 !ELSE ; zpsi = 0. 1635 !ENDIF 1636 ! -- van albada 1 -- 1637 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1638 ! -- smart -- 1639 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1640 ! -- umist -- 1641 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1642 1643 ! high order flux corrected by the limiter 1644 pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 1645 1646 ENDIF 1782 END DO 1647 1783 END DO 1648 1784 END DO 1649 CALL lbc_lnk( pfu_ho, 'U', -1.) ! lateral boundary cond.1785 CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond. 1650 1786 ! 1651 1787 END SUBROUTINE limiter_x … … 1657 1793 !! ** Purpose : compute flux limiter 1658 1794 !!---------------------------------------------------------------------- 1659 REAL(wp) 1660 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pv ! ice i-velocity => u*e21661 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pvc ! ice i-velocity *A => u*e2*a1662 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ) :: pt ! ice tracer1663 REAL(wp), DIMENSION ( jpi,jpj), INTENT(inout) :: pfv_ho ! high order flux1664 REAL(wp), DIMENSION ( jpi,jpj), INTENT(in ), OPTIONAL :: pfv_ups ! upstream flux1795 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1796 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice i-velocity => u*e2 1797 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pvc ! ice i-velocity *A => u*e2*a 1798 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer 1799 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho ! high order flux 1800 REAL(wp), DIMENSION (:,:,:), INTENT(in ), OPTIONAL :: pfv_ups ! upstream flux 1665 1801 ! 1666 1802 REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr 1667 INTEGER :: ji, jj ! dummy loop indices 1668 REAL(wp), DIMENSION (jpi,jpj) :: zslpy ! tracer slopes 1669 !!---------------------------------------------------------------------- 1670 ! 1671 DO jj = 2, jpjm1 1672 DO ji = fs_2, fs_jpim1 ! vector opt. 1673 zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1) 1803 INTEGER :: ji, jj, jl ! dummy loop indices 1804 REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpy ! tracer slopes 1805 !!---------------------------------------------------------------------- 1806 ! 1807 DO jl = 1, jpl 1808 DO jj = 2, jpjm1 1809 DO ji = fs_2, fs_jpim1 ! vector opt. 1810 zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 1811 END DO 1674 1812 END DO 1675 1813 END DO 1676 CALL lbc_lnk( zslpy, 'V', -1.) ! lateral boundary cond. 1677 1678 DO jj = 2, jpjm1 1679 DO ji = fs_2, fs_jpim1 ! vector opt. 1680 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1681 1682 Rjm = zslpy(ji,jj-1) 1683 Rj = zslpy(ji,jj ) 1684 Rjp = zslpy(ji,jj+1) 1685 1686 IF( PRESENT(pfv_ups) ) THEN 1687 1688 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 1689 ELSE ; Rr = Rjp 1814 CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond. 1815 1816 DO jl = 1, jpl 1817 DO jj = 2, jpjm1 1818 DO ji = fs_2, fs_jpim1 ! vector opt. 1819 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1820 1821 Rjm = zslpy(ji,jj-1,jl) 1822 Rj = zslpy(ji,jj ,jl) 1823 Rjp = zslpy(ji,jj+1,jl) 1824 1825 IF( PRESENT(pfv_ups) ) THEN 1826 1827 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 1828 ELSE ; Rr = Rjp 1829 ENDIF 1830 1831 zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1832 IF( Rj > 0. ) THEN 1833 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), & 1834 & MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1835 ELSE 1836 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), & 1837 & MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1838 ENDIF 1839 pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 1840 1841 ELSE 1842 1843 IF( Rj /= 0. ) THEN 1844 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1845 ELSE ; Cr = Rjp / Rj 1846 ENDIF 1847 ELSE 1848 Cr = 0. 1849 !IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 1850 !ELSE ; Cr = Rjp * 1.e20 1851 !ENDIF 1852 ENDIF 1853 1854 ! -- superbee -- 1855 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1856 ! -- van albada 2 -- 1857 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1858 1859 ! -- sweby (with beta=1) -- 1860 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1861 ! -- van Leer -- 1862 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1863 ! -- ospre -- 1864 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1865 ! -- koren -- 1866 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1867 ! -- charm -- 1868 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1869 !ELSE ; zpsi = 0. 1870 !ENDIF 1871 ! -- van albada 1 -- 1872 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1873 ! -- smart -- 1874 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1875 ! -- umist -- 1876 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1877 1878 ! high order flux corrected by the limiter 1879 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1880 1690 1881 ENDIF 1691 1692 zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj) 1693 IF( Rj > 0. ) THEN 1694 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)), & 1695 & MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 1696 ELSE 1697 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)), & 1698 & MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 1699 ENDIF 1700 pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter 1701 1702 ELSE 1703 1704 IF( Rj /= 0. ) THEN 1705 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1706 ELSE ; Cr = Rjp / Rj 1707 ENDIF 1708 ELSE 1709 Cr = 0. 1710 !IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 1711 !ELSE ; Cr = Rjp * 1.e20 1712 !ENDIF 1713 ENDIF 1714 1715 ! -- superbee -- 1716 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1717 ! -- van albada 2 -- 1718 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1719 1720 ! -- sweby (with beta=1) -- 1721 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1722 ! -- van Leer -- 1723 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1724 ! -- ospre -- 1725 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1726 ! -- koren -- 1727 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1728 ! -- charm -- 1729 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1730 !ELSE ; zpsi = 0. 1731 !ENDIF 1732 ! -- van albada 1 -- 1733 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1734 ! -- smart -- 1735 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1736 ! -- umist -- 1737 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1738 1739 ! high order flux corrected by the limiter 1740 pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1741 1742 ENDIF 1882 END DO 1743 1883 END DO 1744 1884 END DO 1745 CALL lbc_lnk( pfv_ho, 'V', -1.) ! lateral boundary cond.1885 CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond. 1746 1886 ! 1747 1887 END SUBROUTINE limiter_y -
NEMO/trunk/src/ICE/icedyn_rdgrft.F90
r10413 r10425 91 91 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 92 92 93 IF( lk_mpp ) CALL mpp_sum (ice_dyn_rdgrft_alloc )94 IF( ice_dyn_rdgrft_alloc /= 0 ) CALL ctl_ warn( 'ice_dyn_rdgrft_alloc: failed to allocate arrays')93 CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc ) 94 IF( ice_dyn_rdgrft_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_alloc: failed to allocate arrays' ) 95 95 ! 96 96 END FUNCTION ice_dyn_rdgrft_alloc … … 269 269 ! 270 270 iter = iter + 1 271 IF( iter > jp_itermax ) CALL ctl_ warn( 'icedyn_rdgrft: non-converging ridging scheme')271 IF( iter > jp_itermax ) CALL ctl_stop( 'STOP', 'icedyn_rdgrft: non-converging ridging scheme' ) 272 272 ! 273 273 END DO … … 793 793 END DO 794 794 END DO 795 CALL lbc_lnk( strength, 'T', 1. )795 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 796 796 ! 797 797 CASE( 2 ) !--- Temporal smoothing … … 814 814 END DO 815 815 END DO 816 CALL lbc_lnk( strength, 'T', 1. )816 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 817 817 ! 818 818 END SELECT -
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
r10415 r10425 146 146 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 147 147 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 148 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence148 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence 149 149 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 150 150 ! ! ocean surface (ssh_m) if ice is not embedded … … 193 193 END DO 194 194 END DO 195 CALL lbc_lnk( zfmask, 'F', 1._wp )195 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp ) 196 196 197 197 ! Lateral boundary conditions on velocity (modify zfmask) … … 220 220 ENDIF 221 221 END DO 222 CALL lbc_lnk( zfmask, 'F', 1._wp )222 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp ) 223 223 224 224 !------------------------------------------------------------------------------! … … 322 322 END DO 323 323 END DO 324 CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )324 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. ) 325 325 ! 326 326 ! !== Landfast ice parameterization ==! … … 343 343 END DO 344 344 END DO 345 CALL lbc_lnk( tau_icebfr(:,:), 'T', 1. )345 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 346 346 ! 347 347 ELSEIF( ln_landfast_home ) THEN !-- Home made … … 370 370 DO jter = 1 , nn_nevp ! loop over jter ! 371 371 ! !----------------------! 372 IF(ln_ctl) THEN ! Convergence test 373 DO jj = 1, jpjm1 374 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 375 zv_ice(:,jj) = v_ice(:,jj) 376 END DO 377 ENDIF 372 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 373 ! 374 !!$ IF(ln_ctl) THEN ! Convergence test 375 !!$ DO jj = 1, jpjm1 376 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 377 !!$ zv_ice(:,jj) = v_ice(:,jj) 378 !!$ END DO 379 !!$ ENDIF 378 380 379 381 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! … … 388 390 END DO 389 391 END DO 390 CALL lbc_lnk( zds, 'F', 1. )392 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. ) 391 393 392 394 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 … … 432 434 END DO 433 435 END DO 434 CALL lbc_lnk( zp_delt, 'T', 1. )436 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 435 437 436 438 DO jj = 1, jpjm1 … … 527 529 END DO 528 530 END DO 529 CALL lbc_lnk( v_ice, 'V', -1. )531 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. ) 530 532 ! 531 533 #if defined key_agrif … … 575 577 END DO 576 578 END DO 577 CALL lbc_lnk( u_ice, 'U', -1. )579 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. ) 578 580 ! 579 581 #if defined key_agrif … … 625 627 END DO 626 628 END DO 627 CALL lbc_lnk( u_ice, 'U', -1. )629 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. ) 628 630 ! 629 631 #if defined key_agrif … … 673 675 END DO 674 676 END DO 675 CALL lbc_lnk( v_ice, 'V', -1. )677 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. ) 676 678 ! 677 679 #if defined key_agrif … … 682 684 ! 683 685 ENDIF 684 685 IF(ln_ctl) THEN ! Convergence test686 DO jj = 2 , jpjm1687 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )688 END DO689 zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )690 IF( lk_mpp ) CALL mpp_max(zresm ) ! max over the global domain691 ENDIF686 687 !!$ IF(ln_ctl) THEN ! Convergence test 688 !!$ DO jj = 2 , jpjm1 689 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 690 !!$ END DO 691 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 692 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 693 !!$ ENDIF 692 694 ! 693 695 ! ! ==================== ! … … 738 740 END DO 739 741 END DO 740 CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )742 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 741 743 742 744 ! --- Store the stress tensor for the next time step --- ! 743 CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )745 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 744 746 pstress1_i (:,:) = zs1 (:,:) 745 747 pstress2_i (:,:) = zs2 (:,:) … … 785 787 END DO 786 788 END DO 787 CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )789 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 788 790 ! 789 791 IF( iom_use('isig1') ) CALL iom_put( "isig1" , zsig1 ) … … 842 844 END DO 843 845 844 CALL lbc_lnk_multi( zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., &846 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., & 845 847 & zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1., & 846 848 & zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1., & 847 849 & zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1. ) 848 850 849 CALL lbc_lnk_multi( zdiag_utau_oi , 'U', -1., zdiag_vtau_oi , 'V', -1., &851 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi , 'U', -1., zdiag_vtau_oi , 'V', -1., & 850 852 & zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., & 851 853 & zdiag_xatrp , 'U', -1., zdiag_ymtrp_ice, 'V', -1., & -
NEMO/trunk/src/ICE/iceforcing.F90
r10069 r10425 83 83 END DO 84 84 END DO 85 CALL lbc_lnk_multi( utau_ice, 'U', -1., vtau_ice, 'V', -1. )85 CALL lbc_lnk_multi( 'iceforcing', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 86 86 ENDIF 87 87 ! -
NEMO/trunk/src/ICE/iceistate.F90
r10413 r10425 414 414 u_ice (:,:) = 0._wp 415 415 v_ice (:,:) = 0._wp 416 ! fields needed for ice_dyn_adv_umx 417 l_split_advumx(1) = .FALSE. 416 418 ! 417 419 !---------------------------------------------- … … 483 485 !!clem: output of initial state should be written here but it is impossible because 484 486 !! the ocean and ice are in the same file 485 !! CALL dia_wri_state( 'output.init' , nit000)487 !! CALL dia_wri_state( 'output.init' ) 486 488 ! 487 489 END SUBROUTINE ice_istate -
NEMO/trunk/src/ICE/icerst.F90
r10069 r10425 69 69 IF(lwp) THEN 70 70 WRITE(numout,*) 71 SELECT CASE ( jprstlib ) 72 CASE DEFAULT 73 WRITE(numout,*) ' open ice restart NetCDF file: ',TRIM(clpath)//clname 74 END SELECT 71 WRITE(numout,*) ' open ice restart NetCDF file: ',TRIM(clpath)//clname 75 72 IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN 76 73 WRITE(numout,*) ' kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp … … 80 77 ENDIF 81 78 ! 82 CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., k iolib = jprstlib, kdlev = jpl )79 CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kdlev = jpl ) 83 80 lrst_ice = .TRUE. 84 81 ENDIF … … 98 95 INTEGER, INTENT(in) :: kt ! number of iteration 99 96 !! 100 INTEGER :: jk ,jl! dummy loop indices97 INTEGER :: jk ! dummy loop indices 101 98 INTEGER :: iter 102 99 CHARACTER(len=25) :: znam … … 118 115 CALL iom_rstput( iter, nitrst, numriw, 'nn_fsbc', REAL( nn_fsbc, wp ) ) ! time-step 119 116 CALL iom_rstput( iter, nitrst, numriw, 'kt_ice' , REAL( iter , wp ) ) ! date 117 CALL iom_delay_rst( 'WRITE', 'ICE', numriw ) ! save only ice delayed global communication variables 120 118 121 119 ! Prognostic variables … … 169 167 !! ** purpose : read restart file 170 168 !!---------------------------------------------------------------------- 171 INTEGER :: jk , jl169 INTEGER :: jk 172 170 LOGICAL :: llok 173 171 INTEGER :: id1 ! local integer 174 INTEGER :: jlibalt = jprstlib175 172 CHARACTER(len=25) :: znam 176 173 CHARACTER(len=2) :: zchar, zchar1 … … 185 182 ENDIF 186 183 187 CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, k iolib = jprstlib, kdlev = jpl )184 CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, kdlev = jpl ) 188 185 189 186 CALL iom_get( numrir, 'nn_fsbc', zfice ) … … 236 233 CALL iom_get( numrir, jpdom_autoglo, 'u_ice', u_ice ) 237 234 CALL iom_get( numrir, jpdom_autoglo, 'v_ice', v_ice ) 235 236 CALL iom_delay_rst( 'READ', 'ICE', numrir ) ! read only ice delayed global communication variables 237 238 238 ! fields needed for Met Office (Jules) coupling 239 239 IF( ln_cpl ) THEN -
NEMO/trunk/src/ICE/icestp.F90
r10415 r10425 241 241 ierr = ierr + ice1D_alloc () ! thermodynamics 242 242 ! 243 IF( lk_mpp ) CALL mpp_sum(ierr )243 CALL mpp_sum( 'icestp', ierr ) 244 244 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays') 245 245 ! -
NEMO/trunk/src/ICE/icethd.F90
r10069 r10425 126 126 END DO 127 127 ENDIF 128 CALL lbc_lnk( zfric, 'T', 1. )128 CALL lbc_lnk( 'icethd', zfric, 'T', 1. ) 129 129 ! 130 130 !--------------------------------------------------------------------! … … 212 212 END DO 213 213 214 IF( lk_mpp ) CALL mpp_ini_ice( npti , numout )215 216 214 IF( npti > 0 ) THEN ! If there is no ice, do nothing. 217 215 ! … … 250 248 ! ! --- & Move to 2D arrays --- ! 251 249 ! 252 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??253 250 ENDIF 254 251 ! -
NEMO/trunk/src/ICE/icethd_do.F90
r10069 r10425 189 189 END DO 190 190 ! 191 CALL lbc_lnk_multi( zvrel, 'T', 1., ht_i_new, 'T', 1. )191 CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1., ht_i_new, 'T', 1. ) 192 192 193 193 ENDIF -
NEMO/trunk/src/ICE/icethd_zdf_bl99.F90
r10069 r10425 83 83 INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation 84 84 INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation 85 ! 85 86 LOGICAL, DIMENSION(jpij) :: l_T_converged ! true when T converges (per grid point) 87 ! 86 88 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 87 89 REAL(wp) :: zg1 = 2._wp ! … … 113 115 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 114 116 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 117 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy 115 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 116 119 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice … … 201 204 ! 202 205 iconv = 0 ! number of iterations 203 zdti_max = 1000._wp ! maximal value of error on all points204 !206 ! 207 l_T_converged(:) = .FALSE. 205 208 ! !============================! 206 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! 209 ! Convergence calculated until all sub-domain grid points have converged 210 ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 211 ! but values are not taken into account (results independant of MPI partitioning) 212 ! 213 DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max ) ! Iterative procedure begins ! 207 214 ! !============================! 208 215 iconv = iconv + 1 … … 217 224 ! 218 225 DO ji = 1, npti 219 ztcond_i (ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )220 ztcond_i (ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )226 ztcond_i_cp(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 227 ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 221 228 END DO 222 229 DO jk = 1, nlay_i-1 223 230 DO ji = 1, npti 224 ztcond_i (ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &225 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )231 ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 232 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 226 233 END DO 227 234 END DO … … 230 237 ! 231 238 DO ji = 1, npti 232 ztcond_i (ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &233 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 )234 ztcond_i (ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &235 & - 0.011_wp * ( t_bo_1d(ji) - rt0 )239 ztcond_i_cp(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 240 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 241 ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 242 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 236 243 END DO 237 244 DO jk = 1, nlay_i-1 238 245 DO ji = 1, npti 239 ztcond_i (ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &240 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) &241 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )246 ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 247 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) & 248 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 242 249 END DO 243 250 END DO 244 251 ! 245 252 ENDIF 246 ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) ) 253 254 ! Variable used after iterations 255 ! Value must be frozen after convergence for MPP independance reason 256 DO ji = 1, npti 257 IF ( .NOT. l_T_converged(ji) ) & 258 ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) ) 259 END DO 247 260 ! 248 261 !--- G(he) : enhancement of thermal conductivity in mono-category case … … 270 283 !----------------- 271 284 !--- Snow 285 ! Variable used after iterations 286 ! Value must be frozen after convergence for MPP independance reason 272 287 DO jk = 0, nlay_s-1 273 288 DO ji = 1, npti 274 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 289 IF ( .NOT. l_T_converged(ji) ) & 290 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 275 291 END DO 276 292 END DO 277 293 DO ji = 1, npti ! Snow-ice interface 278 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 279 IF( zfac > epsi10 ) THEN 280 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 281 ELSE 282 zkappa_s(ji,nlay_s) = 0._wp 294 IF ( .NOT. l_T_converged(ji) ) THEN 295 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 296 IF( zfac > epsi10 ) THEN 297 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 298 ELSE 299 zkappa_s(ji,nlay_s) = 0._wp 300 ENDIF 283 301 ENDIF 284 302 END DO 285 303 286 304 !--- Ice 305 ! Variable used after iterations 306 ! Value must be frozen after convergence for MPP independance reason 287 307 DO jk = 0, nlay_i 288 308 DO ji = 1, npti 289 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 309 IF ( .NOT. l_T_converged(ji) ) & 310 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 290 311 END DO 291 312 END DO 292 313 DO ji = 1, npti ! Snow-ice interface 293 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 314 IF ( .NOT. l_T_converged(ji) ) & 315 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 294 316 END DO 295 317 ! … … 326 348 ! update of the non solar flux according to the update in T_su 327 349 DO ji = 1, npti 328 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 350 ! Variable used after iterations 351 ! Value must be frozen after convergence for MPP independance reason 352 IF ( .NOT. l_T_converged(ji) ) & 353 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 329 354 END DO 330 355 … … 496 521 ! ice temperatures 497 522 DO ji = 1, npti 498 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 523 ! Variable used after iterations 524 ! Value must be frozen after convergence for MPP independance reason 525 IF ( .NOT. l_T_converged(ji) ) & 526 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 499 527 END DO 500 528 … … 502 530 DO ji = 1, npti 503 531 jk = jm - nlay_s - 1 504 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 532 IF ( .NOT. l_T_converged(ji) ) & 533 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 505 534 END DO 506 535 END DO 507 536 508 537 DO ji = 1, npti 509 ! snow temperatures 510 IF( h_s_1d(ji) > 0._wp ) THEN 511 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 512 ENDIF 513 ! surface temperature 514 ztsub(ji) = t_su_1d(ji) 515 IF( t_su_1d(ji) < rt0 ) THEN 516 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 517 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 538 ! Variables used after iterations 539 ! Value must be frozen after convergence for MPP independance reason 540 IF ( .NOT. l_T_converged(ji) ) THEN 541 ! snow temperatures 542 IF( h_s_1d(ji) > 0._wp ) THEN 543 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 544 ENDIF 545 ! surface temperature 546 ztsub(ji) = t_su_1d(ji) 547 IF( t_su_1d(ji) < rt0 ) THEN 548 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 549 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 550 ENDIF 518 551 ENDIF 519 552 END DO … … 524 557 ! check that nowhere it has started to melt 525 558 ! zdti_max is a measure of error, it has to be under zdti_bnd 526 zdti_max = 0._wp 527 DO ji = 1, npti 528 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 529 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 530 END DO 531 532 DO jk = 1, nlay_s 533 DO ji = 1, npti 534 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 535 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 536 END DO 537 END DO 538 539 DO jk = 1, nlay_i 540 DO ji = 1, npti 541 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 542 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 543 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 544 END DO 545 END DO 546 547 ! Compute spatial maximum over all errors 548 ! note that this could be optimized substantially by iterating only the non-converging points 549 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 550 ! 559 560 DO ji = 1, npti 561 562 zdti_max = 0._wp 563 564 IF ( .NOT. l_T_converged(ji) ) THEN 565 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 566 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 567 568 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 569 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 570 571 DO jk = 1, nlay_i 572 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 573 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 574 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 575 END DO 576 577 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 578 579 ENDIF 580 581 END DO 582 551 583 !----------------------------------------! 552 584 ! ! … … 670 702 671 703 ! ice temperatures 672 DO ji = 1, npti 673 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 704 DO ji = 1, npti 705 ! Variable used after iterations 706 ! Value must be frozen after convergence for MPP independance reason 707 IF ( .NOT. l_T_converged(ji) ) & 708 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 674 709 END DO 675 710 676 711 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 677 712 DO ji = 1, npti 678 jk = jm - nlay_s - 1 679 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 713 IF ( .NOT. l_T_converged(ji) ) THEN 714 jk = jm - nlay_s - 1 715 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 716 ENDIF 680 717 END DO 681 718 END DO … … 683 720 ! snow temperatures 684 721 DO ji = 1, npti 685 IF( h_s_1d(ji) > 0._wp ) THEN 686 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 722 ! Variable used after iterations 723 ! Value must be frozen after convergence for MPP independance reason 724 IF ( .NOT. l_T_converged(ji) ) THEN 725 IF( h_s_1d(ji) > 0._wp ) THEN 726 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 727 ENDIF 687 728 ENDIF 688 729 END DO … … 693 734 ! check that nowhere it has started to melt 694 735 ! zdti_max is a measure of error, it has to be under zdti_bnd 695 zdti_max = 0._wp 696 697 DO jk = 1, nlay_s 698 DO ji = 1, npti 699 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 700 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 701 END DO 702 END DO 703 704 DO jk = 1, nlay_i 705 DO ji = 1, npti 706 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 707 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 708 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 709 END DO 710 END DO 711 712 ! Compute spatial maximum over all errors 713 ! note that this could be optimized substantially by iterating only the non-converging points 714 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 736 737 DO ji = 1, npti 738 739 zdti_max = 0._wp 740 741 IF ( .NOT. l_T_converged(ji) ) THEN 742 ! t_s 743 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 744 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 745 ! t_i 746 DO jk = 0, nlay_i 747 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 748 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 749 zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 750 END DO 751 752 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 753 754 ENDIF 755 756 END DO 715 757 716 758 ENDIF ! k_jules -
NEMO/trunk/src/ICE/iceupdate.F90
r10069 r10425 59 59 ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc ) 60 60 ! 61 IF( lk_mpp ) CALL mpp_sum(ice_update_alloc )62 IF( ice_update_alloc /= 0 ) CALL ctl_ warn('ice_update_alloc: failed to allocate arrays')61 CALL mpp_sum( 'iceupdate', ice_update_alloc ) 62 IF( ice_update_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_update_alloc: failed to allocate arrays' ) 63 63 ! 64 64 END FUNCTION ice_update_alloc … … 350 350 END DO 351 351 END DO 352 CALL lbc_lnk_multi( taum, 'T', 1., tmod_io, 'T', 1. )352 CALL lbc_lnk_multi( 'iceupdate', taum, 'T', 1., tmod_io, 'T', 1. ) 353 353 ! 354 354 utau_oce(:,:) = utau(:,:) !* save the air-ocean stresses at ice time-step … … 374 374 END DO 375 375 END DO 376 CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. ) ! lateral boundary condition376 CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1., vtau, 'V', -1. ) ! lateral boundary condition 377 377 ! 378 378 IF( ln_timing ) CALL timing_stop('ice_update_tau') -
NEMO/trunk/src/ICE/icewri.F90
r10413 r10425 141 141 END DO 142 142 END DO 143 CALL lbc_lnk( z2d, 'T', 1. )143 CALL lbc_lnk( 'icewri', z2d, 'T', 1. ) 144 144 IF( iom_use('icevel') ) CALL iom_put( "icevel" , z2d ) 145 145 … … 191 191 ELSEWHERE ; zmsk00(:,:) = 0. 192 192 END WHERE 193 zdiag_area_nh = glob_sum( at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )194 zdiag_volu_nh = glob_sum( vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )193 zdiag_area_nh = glob_sum( 'icewri', at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 194 zdiag_volu_nh = glob_sum( 'icewri', vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 195 195 ! 196 196 WHERE( ff_t > 0._wp .AND. at_i > 0.15 ) ; zmsk00(:,:) = 1.0e-12 197 197 ELSEWHERE ; zmsk00(:,:) = 0. 198 198 END WHERE 199 zdiag_extt_nh = glob_sum( zmsk00(:,:) * e1e2t(:,:) )199 zdiag_extt_nh = glob_sum( 'icewri', zmsk00(:,:) * e1e2t(:,:) ) 200 200 ! 201 201 IF( iom_use('NH_icearea') ) CALL iom_put( "NH_icearea" , zdiag_area_nh ) … … 210 210 ELSEWHERE ; zmsk00(:,:) = 0. 211 211 END WHERE 212 zdiag_area_sh = glob_sum( at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )213 zdiag_volu_sh = glob_sum( vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )212 zdiag_area_sh = glob_sum( 'icewri', at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 213 zdiag_volu_sh = glob_sum( 'icewri', vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 214 214 ! 215 215 WHERE( ff_t < 0._wp .AND. at_i > 0.15 ); zmsk00(:,:) = 1.0e-12 216 216 ELSEWHERE ; zmsk00(:,:) = 0. 217 217 END WHERE 218 zdiag_extt_sh = glob_sum( zmsk00(:,:) * e1e2t(:,:) )218 zdiag_extt_sh = glob_sum( 'icewri', zmsk00(:,:) * e1e2t(:,:) ) 219 219 ! 220 220 IF( iom_use('SH_icearea') ) CALL iom_put( "SH_icearea", zdiag_area_sh ) … … 234 234 235 235 236 SUBROUTINE ice_wri_state( k t, kid, kh_i)236 SUBROUTINE ice_wri_state( kid ) 237 237 !!--------------------------------------------------------------------- 238 238 !! *** ROUTINE ice_wri_state *** … … 245 245 !! History : 4.0 ! 2013-06 (C. Rousset) 246 246 !!---------------------------------------------------------------------- 247 INTEGER, INTENT( in ) :: kt ! ocean time-step index 248 INTEGER, INTENT( in ) :: kid , kh_i 249 INTEGER :: nz_i, jl 250 REAL(wp), DIMENSION(jpl) :: jcat 247 INTEGER, INTENT( in ) :: kid 251 248 !!---------------------------------------------------------------------- 252 249 ! 253 DO jl = 1, jpl 254 jcat(jl) = REAL(jl) 255 END DO 256 257 CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up") 258 259 CALL histdef( kid, "sithic", "Ice thickness" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 260 CALL histdef( kid, "siconc", "Ice concentration" , "%" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 261 CALL histdef( kid, "sitemp", "Ice temperature" , "C" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 262 CALL histdef( kid, "sivelu", "i-Ice speed " , "m/s" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 263 CALL histdef( kid, "sivelv", "j-Ice speed " , "m/s" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 264 CALL histdef( kid, "sistru", "i-Wind stress over ice" , "Pa" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 265 CALL histdef( kid, "sistrv", "j-Wind stress over ice" , "Pa" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 266 CALL histdef( kid, "sisflx", "Solar flx over ocean" , "W/m2" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 267 CALL histdef( kid, "sinflx", "NonSolar flx over ocean", "W/m2" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 268 CALL histdef( kid, "snwpre", "Snow precipitation" , "kg/m2/s", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 269 CALL histdef( kid, "sisali", "Ice salinity" , "PSU" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 270 CALL histdef( kid, "sivolu", "Ice volume" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 271 CALL histdef( kid, "sidive", "Ice divergence" , "10-8s-1", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 272 CALL histdef( kid, "si_amp", "Melt pond fraction" , "%" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 273 CALL histdef( kid, "si_vmp", "Melt pond volume" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 274 ! 275 CALL histdef( kid, "sithicat", "Ice thickness" , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 276 CALL histdef( kid, "siconcat", "Ice concentration" , "%" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 277 CALL histdef( kid, "sisalcat", "Ice salinity" , "" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 278 CALL histdef( kid, "snthicat", "Snw thickness" , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 279 280 CALL histend( kid, snc4set ) ! end of the file definition 281 282 CALL histwrite( kid, "sithic", kt, hm_i , jpi*jpj, (/1/) ) 283 CALL histwrite( kid, "siconc", kt, at_i , jpi*jpj, (/1/) ) 284 CALL histwrite( kid, "sitemp", kt, tm_i - rt0 , jpi*jpj, (/1/) ) 285 CALL histwrite( kid, "sivelu", kt, u_ice , jpi*jpj, (/1/) ) 286 CALL histwrite( kid, "sivelv", kt, v_ice , jpi*jpj, (/1/) ) 287 CALL histwrite( kid, "sistru", kt, utau_ice , jpi*jpj, (/1/) ) 288 CALL histwrite( kid, "sistrv", kt, vtau_ice , jpi*jpj, (/1/) ) 289 CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) ) 290 CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) ) 291 CALL histwrite( kid, "snwpre", kt, sprecip , jpi*jpj, (/1/) ) 292 CALL histwrite( kid, "sisali", kt, sm_i , jpi*jpj, (/1/) ) 293 CALL histwrite( kid, "sivolu", kt, vt_i , jpi*jpj, (/1/) ) 294 CALL histwrite( kid, "sidive", kt, divu_i*1.0e8 , jpi*jpj, (/1/) ) 295 CALL histwrite( kid, "si_amp", kt, at_ip , jpi*jpj, (/1/) ) 296 CALL histwrite( kid, "si_vmp", kt, vt_ip , jpi*jpj, (/1/) ) 297 ! 298 CALL histwrite( kid, "sithicat", kt, h_i , jpi*jpj*jpl, (/1/) ) 299 CALL histwrite( kid, "siconcat", kt, a_i , jpi*jpj*jpl, (/1/) ) 300 CALL histwrite( kid, "sisalcat", kt, s_i , jpi*jpj*jpl, (/1/) ) 301 CALL histwrite( kid, "snthicat", kt, h_s , jpi*jpj*jpl, (/1/) ) 302 303 !! The file is closed in dia_wri_state (ocean routine) 304 !! CALL histclo( kid ) 305 ! 250 !! The file is open in dia_wri_state (ocean routine) 251 252 CALL iom_rstput( 0, 0, kid, 'sithic', hm_i ) ! Ice thickness 253 CALL iom_rstput( 0, 0, kid, 'siconc', at_i ) ! Ice concentration 254 CALL iom_rstput( 0, 0, kid, 'sitemp', tm_i - rt0 ) ! Ice temperature 255 CALL iom_rstput( 0, 0, kid, 'sivelu', u_ice ) ! i-Ice speed 256 CALL iom_rstput( 0, 0, kid, 'sivelv', v_ice ) ! j-Ice speed 257 CALL iom_rstput( 0, 0, kid, 'sistru', utau_ice ) ! i-Wind stress over ice 258 CALL iom_rstput( 0, 0, kid, 'sistrv', vtau_ice ) ! i-Wind stress over ice 259 CALL iom_rstput( 0, 0, kid, 'sisflx', qsr ) ! Solar flx over ocean 260 CALL iom_rstput( 0, 0, kid, 'sinflx', qns ) ! NonSolar flx over ocean 261 CALL iom_rstput( 0, 0, kid, 'snwpre', sprecip ) ! Snow precipitation 262 CALL iom_rstput( 0, 0, kid, 'sisali', sm_i ) ! Ice salinity 263 CALL iom_rstput( 0, 0, kid, 'sivolu', vt_i ) ! Ice volume 264 CALL iom_rstput( 0, 0, kid, 'sidive', divu_i*1.0e8 ) ! Ice divergence 265 CALL iom_rstput( 0, 0, kid, 'si_amp', at_ip ) ! Melt pond fraction 266 CALL iom_rstput( 0, 0, kid, 'si_vmp', vt_ip ) ! Melt pond volume 267 CALL iom_rstput( 0, 0, kid, 'sithicat', h_i ) ! Ice thickness 268 CALL iom_rstput( 0, 0, kid, 'siconcat', a_i ) ! Ice concentration 269 CALL iom_rstput( 0, 0, kid, 'sisalcat', s_i ) ! Ice salinity 270 CALL iom_rstput( 0, 0, kid, 'snthicat', h_s ) ! Snw thickness 271 306 272 END SUBROUTINE ice_wri_state 307 273
Note: See TracChangeset
for help on using the changeset viewer.
