Changeset 11612 for NEMO/trunk
- Timestamp:
- 2019-09-29T10:55:44+02:00 (5 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 5 edited
-
cfgs/SHARED/namelist_ice_ref (modified) (1 diff)
-
cfgs/SPITZ12/EXPREF/namelist_ice_cfg (modified) (1 diff)
-
src/ICE/icectl.F90 (modified) (1 diff)
-
src/ICE/icedyn_adv_pra.F90 (modified) (15 diffs)
-
src/ICE/icedyn_adv_umx.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/cfgs/SHARED/namelist_ice_ref
r11601 r11612 103 103 &namdyn_adv ! Ice advection 104 104 !------------------------------------------------------------------------------ 105 ln_adv_Pra = . false. ! Advection scheme (Prather)106 ln_adv_UMx = . true. ! Advection scheme (Ultimate-Macho)105 ln_adv_Pra = .true. ! Advection scheme (Prather) 106 ln_adv_UMx = .false. ! Advection scheme (Ultimate-Macho) 107 107 nn_UMx = 5 ! order of the scheme for UMx (1-5 ; 20=centered 2nd order) 108 108 / -
NEMO/trunk/cfgs/SPITZ12/EXPREF/namelist_ice_cfg
r11548 r11612 50 50 &namdyn_adv ! Ice advection 51 51 !------------------------------------------------------------------------------ 52 ln_adv_Pra = .false. ! Advection scheme (Prather) 53 ln_adv_UMx = .true. ! Advection scheme (Ultimate-Macho) 54 nn_UMx = 5 ! order of the scheme for UMx (1-5 ; 20=centered 2nd order) 52 55 / 53 56 !------------------------------------------------------------------------------ -
NEMO/trunk/src/ICE/icectl.F90
r11601 r11612 158 158 ! only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 159 159 ! so the formulation for conservation is different (and not coded) 160 ! it does not mean UM is not conservative (it is checked with above prints) 161 IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &162 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice160 ! it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 161 !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 162 ! & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 163 163 ENDIF 164 164 ! -
NEMO/trunk/src/ICE/icedyn_adv_pra.F90
r10425 r11612 19 19 USE ice ! sea-ice variables 20 20 USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call 21 USE icevar ! sea-ice: operations 21 22 ! 22 23 USE in_out_manager ! I/O manager … … 25 26 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 26 27 USE lbclnk ! lateral boundary conditions (or mpp links) 27 USE prtctl ! Print control28 28 29 29 IMPLICIT NONE … … 39 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsal, sysal, sxxsal, syysal, sxysal ! ice salinity 40 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxage, syage, sxxage, syyage, sxyage ! ice age 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ):: sxopw, syopw, sxxopw, syyopw, sxyopw ! open water in sea ice41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxopw, syopw, sxxopw, syyopw, sxyopw ! open water in sea ice 42 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxc0 , syc0 , sxxc0 , syyc0 , sxyc0 ! snow layers heat content 43 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxe , sye , sxxe , syye , sxye ! ice layers heat content … … 82 82 ! 83 83 INTEGER :: jk, jl, jt ! dummy loop indices 84 INTEGER :: initad ! number of sub-timestep for the advection 85 REAL(wp) :: zcfl , zusnit ! - - 86 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zarea 87 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0opw 88 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0smi, z0oi 89 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ap , z0vp 90 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z0es 91 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z0ei 84 INTEGER :: icycle ! number of sub-timestep for the advection 85 REAL(wp) :: zdt ! - - 86 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 87 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 88 REAL(wp), DIMENSION(jpi,jpj,1) :: z0opw 89 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 90 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp 91 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 92 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei 92 93 !!---------------------------------------------------------------------- 93 94 ! 94 95 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 95 96 ! 96 ALLOCATE( zarea(jpi,jpj) , z0opw(jpi,jpj, 1 ) , z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , & 97 & z0ai(jpi,jpj,jpl) , z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap (jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) , & 98 & z0es (jpi,jpj,nlay_s,jpl), z0ei(jpi,jpj,nlay_i,jpl) ) 99 ! 100 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 101 zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 102 zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 103 CALL mpp_max( 'icedyn_adv_pra', zcfl ) 97 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! 98 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 99 ! this should not affect too much the stability 100 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 101 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 104 102 105 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 106 ELSE ; initad = 1 ; zusnit = 1.0_wp 103 ! non-blocking global communication send zcflnow and receive zcflprv 104 CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 105 106 IF( zcflprv(1) > .5 ) THEN ; icycle = 2 107 ELSE ; icycle = 1 107 108 ENDIF 109 zdt = rdt_ice / REAL(icycle) 108 110 109 zarea(:,:) = e1e2t(:,:)110 111 !------------------------- 111 112 ! transported fields … … 113 114 z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:) ! Open water area 114 115 DO jl = 1, jpl 115 z0snw(:,:,jl) = pv_s (:,:, jl) * e1e2t(:,:) ! Snow volume 116 z0ice(:,:,jl) = pv_i (:,:, jl) * e1e2t(:,:) ! Ice volume 117 z0ai (:,:,jl) = pa_i (:,:, jl) * e1e2t(:,:) ! Ice area 118 z0smi(:,:,jl) = psv_i(:,:, jl) * e1e2t(:,:) ! Salt content 119 z0oi (:,:,jl) = poa_i(:,:, jl) * e1e2t(:,:) ! Age content 116 zarea(:,:,jl) = e1e2t(:,:) 117 z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume 118 z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume 119 z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area 120 z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content 121 z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content 120 122 DO jk = 1, nlay_s 121 123 z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content … … 133 135 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 134 136 ! !--------------------------------------------! 135 DO jt = 1, initad 136 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 137 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 138 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 139 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 140 DO jl = 1, jpl 141 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 142 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 143 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 144 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 145 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 146 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 147 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 148 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 149 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 150 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 151 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 152 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 153 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 154 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 155 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 156 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 157 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 158 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 159 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 160 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 161 DO jk = 1, nlay_s !--- snow heat contents --- 162 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), & 163 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 164 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), & 165 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 166 END DO 167 DO jk = 1, nlay_i !--- ice heat contents --- 168 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), & 169 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 170 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), & 171 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 172 END DO 173 IF ( ln_pnd_H12 ) THEN 174 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction -- 175 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 176 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & 177 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 178 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume -- 179 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 180 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & 181 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 182 ENDIF 183 END DO 137 DO jt = 1, icycle 138 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water 139 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) 140 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 141 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 142 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 143 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 144 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 145 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 146 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 147 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 148 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 149 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 150 ! 151 DO jk = 1, nlay_s !--- snow heat content 152 CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 153 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 154 CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 155 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 156 END DO 157 DO jk = 1, nlay_i !--- ice heat content 158 CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 159 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 160 CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 161 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 162 END DO 163 ! 164 IF ( ln_pnd_H12 ) THEN 165 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 166 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 167 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 168 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 169 ENDIF 184 170 END DO 185 171 ! !--------------------------------------------! 186 172 ELSE !== even ice time step: adv_y then adv_x ==! 187 173 ! !--------------------------------------------! 188 DO jt = 1, initad 189 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 190 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 191 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 192 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 193 DO jl = 1, jpl 194 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 195 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 196 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 197 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 198 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 199 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 200 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 201 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 202 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 203 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 204 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 205 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 206 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 207 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 208 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 209 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 210 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 211 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 212 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 213 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 214 DO jk = 1, nlay_s !--- snow heat contents --- 215 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), & 216 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 217 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), & 218 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 219 END DO 220 DO jk = 1, nlay_i !--- ice heat contents --- 221 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), & 222 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 223 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), & 224 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 225 END DO 226 IF ( ln_pnd_H12 ) THEN 227 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction --- 228 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 229 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & 230 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 231 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume --- 232 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 233 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & 234 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 235 ENDIF 236 END DO 174 DO jt = 1, icycle 175 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water 176 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) 177 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 178 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 179 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 180 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 181 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 182 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 183 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 184 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 185 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 186 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 187 DO jk = 1, nlay_s !--- snow heat content 188 CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 189 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 190 CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 191 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 192 END DO 193 DO jk = 1, nlay_i !--- ice heat content 194 CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 195 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 196 CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 197 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 198 END DO 199 IF ( ln_pnd_H12 ) THEN 200 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 201 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 202 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 203 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 204 ENDIF 237 205 END DO 238 206 ENDIF … … 243 211 pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1) 244 212 DO jl = 1, jpl 245 pv_i (:,:, 246 pv_s (:,:, 247 psv_i(:,:, 248 poa_i(:,:, 249 pa_i (:,:, 213 pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 214 pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 215 psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 216 poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 217 pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 250 218 DO jk = 1, nlay_s 251 219 pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) … … 255 223 END DO 256 224 IF ( ln_pnd_H12 ) THEN 257 pa_ip (:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)258 pv_ip (:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)225 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 226 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 259 227 ENDIF 260 228 END DO 261 229 ! 262 DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei ) 230 ! --- Ensure non-negative fields --- ! 231 ! Remove negative values (conservation is ensured) 232 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 233 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 263 234 ! 264 235 IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file … … 267 238 268 239 269 SUBROUTINE adv_x( pd f, put , pcrh, psm , ps0 , &240 SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , & 270 241 & psx, psxx, psy , psyy, psxy ) 271 242 !!---------------------------------------------------------------------- … … 275 246 !! variable on x axis 276 247 !!---------------------------------------------------------------------- 277 REAL(wp) , INTENT(in ) :: pdf ! reduction factor forthe time step278 REAL(wp) 279 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s]280 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psm ! area281 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: ps0 ! field to be advected282 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments283 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments248 REAL(wp) , INTENT(in ) :: pdt ! the time step 249 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 250 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 251 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 252 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 253 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 254 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 284 255 !! 285 INTEGER :: ji, jj 286 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! local scalars256 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 257 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 287 258 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 288 259 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 291 262 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 292 263 !----------------------------------------------------------------------- 293 294 ! Limitation of moments. 295 296 zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 297 298 DO jj = 1, jpj 299 DO ji = 1, jpi 300 zslpmax = MAX( 0._wp, ps0(ji,jj) ) 301 zs1max = 1.5 * zslpmax 302 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) ) 303 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 304 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) ) ) 305 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 306 307 ps0 (ji,jj) = zslpmax 308 psx (ji,jj) = zs1new * rswitch 309 psxx(ji,jj) = zs2new * rswitch 310 psy (ji,jj) = psy (ji,jj) * rswitch 311 psyy(ji,jj) = psyy(ji,jj) * rswitch 312 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 313 END DO 264 ! 265 jcat = SIZE( ps0 , 3 ) ! size of input arrays 266 ! 267 DO jl = 1, jcat ! loop on categories 268 ! 269 ! Limitation of moments. 270 DO jj = 2, jpjm1 271 DO ji = 1, jpi 272 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 273 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 274 ! 275 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 276 zs1max = 1.5 * zslpmax 277 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 278 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 279 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 280 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 281 282 ps0 (ji,jj,jl) = zslpmax 283 psx (ji,jj,jl) = zs1new * rswitch 284 psxx(ji,jj,jl) = zs2new * rswitch 285 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 286 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 287 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 288 END DO 289 END DO 290 291 ! Calculate fluxes and moments between boxes i<-->i+1 292 DO jj = 2, jpjm1 ! Flux from i to i+1 WHEN u GT 0 293 DO ji = 1, jpi 294 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 295 zalf = MAX( 0._wp, put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji,jj,jl) 296 zalfq = zalf * zalf 297 zalf1 = 1.0 - zalf 298 zalf1q = zalf1 * zalf1 299 ! 300 zfm (ji,jj) = zalf * psm (ji,jj,jl) 301 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 302 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 303 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 304 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 305 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 306 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 307 308 ! Readjust moments remaining in the box. 309 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 310 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 311 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 312 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 313 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 314 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 315 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 316 END DO 317 END DO 318 319 DO jj = 2, jpjm1 ! Flux from i+1 to i when u LT 0. 320 DO ji = 1, fs_jpim1 321 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji+1,jj,jl) 322 zalg (ji,jj) = zalf 323 zalfq = zalf * zalf 324 zalf1 = 1.0 - zalf 325 zalg1 (ji,jj) = zalf1 326 zalf1q = zalf1 * zalf1 327 zalg1q(ji,jj) = zalf1q 328 ! 329 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 330 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 331 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 332 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 333 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 334 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 335 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 336 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 337 END DO 338 END DO 339 340 DO jj = 2, jpjm1 ! Readjust moments remaining in the box. 341 DO ji = fs_2, fs_jpim1 342 zbt = zbet(ji-1,jj) 343 zbt1 = 1.0 - zbet(ji-1,jj) 344 ! 345 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 346 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 347 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 348 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 349 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 350 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 351 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 352 END DO 353 END DO 354 355 ! Put the temporary moments into appropriate neighboring boxes. 356 DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0. 357 DO ji = fs_2, fs_jpim1 358 zbt = zbet(ji-1,jj) 359 zbt1 = 1.0 - zbet(ji-1,jj) 360 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 361 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 362 zalf1 = 1.0 - zalf 363 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 364 ! 365 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 366 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 367 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 368 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 369 & + zbt1 * psxx(ji,jj,jl) 370 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 371 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 372 & + zbt1 * psxy(ji,jj,jl) 373 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 374 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 375 END DO 376 END DO 377 378 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0. 379 DO ji = fs_2, fs_jpim1 380 zbt = zbet(ji,jj) 381 zbt1 = 1.0 - zbet(ji,jj) 382 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 383 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 384 zalf1 = 1.0 - zalf 385 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 386 ! 387 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 388 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 389 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 390 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 391 & + ( zalf1 - zalf ) * ztemp ) ) 392 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 393 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 394 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 395 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 396 END DO 397 END DO 398 314 399 END DO 315 400 316 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)317 psm (:,:) = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )318 319 ! Calculate fluxes and moments between boxes i<-->i+1320 DO jj = 1, jpj ! Flux from i to i+1 WHEN u GT 0321 DO ji = 1, jpi322 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )323 zalf = MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)324 zalfq = zalf * zalf325 zalf1 = 1.0 - zalf326 zalf1q = zalf1 * zalf1327 !328 zfm (ji,jj) = zalf * psm (ji,jj)329 zf0 (ji,jj) = zalf * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) ) )330 zfx (ji,jj) = zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )331 zfxx(ji,jj) = zalf * psxx(ji,jj) * zalfq332 zfy (ji,jj) = zalf * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )333 zfxy(ji,jj) = zalfq * psxy(ji,jj)334 zfyy(ji,jj) = zalf * psyy(ji,jj)335 336 ! Readjust moments remaining in the box.337 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj)338 ps0 (ji,jj) = ps0 (ji,jj) - zf0(ji,jj)339 psx (ji,jj) = zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )340 psxx(ji,jj) = zalf1 * zalf1q * psxx(ji,jj)341 psy (ji,jj) = psy (ji,jj) - zfy(ji,jj)342 psyy(ji,jj) = psyy(ji,jj) - zfyy(ji,jj)343 psxy(ji,jj) = zalf1q * psxy(ji,jj)344 END DO345 END DO346 347 DO jj = 1, jpjm1 ! Flux from i+1 to i when u LT 0.348 DO ji = 1, fs_jpim1349 zalf = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj)350 zalg (ji,jj) = zalf351 zalfq = zalf * zalf352 zalf1 = 1.0 - zalf353 zalg1 (ji,jj) = zalf1354 zalf1q = zalf1 * zalf1355 zalg1q(ji,jj) = zalf1q356 !357 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj)358 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )359 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )360 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj) * zalfq361 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )362 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj)363 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj)364 END DO365 END DO366 367 DO jj = 2, jpjm1 ! Readjust moments remaining in the box.368 DO ji = fs_2, fs_jpim1369 zbt = zbet(ji-1,jj)370 zbt1 = 1.0 - zbet(ji-1,jj)371 !372 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )373 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )374 psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )375 psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)376 psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )377 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )378 psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)379 END DO380 END DO381 382 ! Put the temporary moments into appropriate neighboring boxes.383 DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0.384 DO ji = fs_2, fs_jpim1385 zbt = zbet(ji-1,jj)386 zbt1 = 1.0 - zbet(ji-1,jj)387 psm(ji,jj) = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)388 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj)389 zalf1 = 1.0 - zalf390 ztemp = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)391 !392 ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)393 psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)394 psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj) &395 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &396 & + zbt1 * psxx(ji,jj)397 psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj) &398 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj) ) ) &399 & + zbt1 * psxy(ji,jj)400 psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)401 psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)402 END DO403 END DO404 405 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0.406 DO ji = fs_2, fs_jpim1407 zbt = zbet(ji,jj)408 zbt1 = 1.0 - zbet(ji,jj)409 psm(ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )410 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj)411 zalf1 = 1.0 - zalf412 ztemp = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)413 !414 ps0(ji,jj) = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )415 psx(ji,jj) = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )416 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj) &417 & + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) &418 & + ( zalf1 - zalf ) * ztemp ) )419 psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) &420 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) ) )421 psy(ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )422 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )423 END DO424 END DO425 426 401 !-- Lateral boundary conditions 427 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T', 1., ps0 , 'T', 1. & 428 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 429 & , psxx, 'T', 1., psyy, 'T', 1. & 430 & , psxy, 'T', 1. ) 431 432 IF(ln_ctl) THEN 433 CALL prt_ctl(tab2d_1=psm , clinfo1=' adv_x: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ') 434 CALL prt_ctl(tab2d_1=psx , clinfo1=' adv_x: psx :', tab2d_2=psxx, clinfo2=' psxx : ') 435 CALL prt_ctl(tab2d_1=psy , clinfo1=' adv_x: psy :', tab2d_2=psyy, clinfo2=' psyy : ') 436 CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :') 437 ENDIF 402 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 403 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 404 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 438 405 ! 439 406 END SUBROUTINE adv_x 440 407 441 408 442 SUBROUTINE adv_y( pd f, pvt , pcrh, psm , ps0 , &409 SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , & 443 410 & psx, psxx, psy , psyy, psxy ) 444 411 !!--------------------------------------------------------------------- … … 448 415 !! variable on y axis 449 416 !!--------------------------------------------------------------------- 450 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for thetime step451 REAL(wp) 452 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s]453 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psm ! area454 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: ps0 ! field to be advected455 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments456 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments417 REAL(wp) , INTENT(in ) :: pdt ! time step 418 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 419 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 420 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 421 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 422 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 423 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 457 424 !! 458 INTEGER :: ji, jj 459 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! temporary scalars425 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 426 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 460 427 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 461 428 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 464 431 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 465 432 !--------------------------------------------------------------------- 466 467 ! Limitation of moments. 468 469 zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 470 471 DO jj = 1, jpj 472 DO ji = 1, jpi 473 zslpmax = MAX( 0._wp, ps0(ji,jj) ) 474 zs1max = 1.5 * zslpmax 475 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 476 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 477 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) ) ) 478 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 479 ! 480 ps0 (ji,jj) = zslpmax 481 psx (ji,jj) = psx (ji,jj) * rswitch 482 psxx(ji,jj) = psxx(ji,jj) * rswitch 483 psy (ji,jj) = zs1new * rswitch 484 psyy(ji,jj) = zs2new * rswitch 485 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 486 END DO 433 ! 434 jcat = SIZE( ps0 , 3 ) ! size of input arrays 435 ! 436 DO jl = 1, jcat ! loop on categories 437 ! 438 ! Limitation of moments. 439 DO jj = 1, jpj 440 DO ji = fs_2, fs_jpim1 441 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 442 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 443 ! 444 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 445 zs1max = 1.5 * zslpmax 446 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 447 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 448 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 449 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 450 ! 451 ps0 (ji,jj,jl) = zslpmax 452 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 453 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 454 psy (ji,jj,jl) = zs1new * rswitch 455 psyy(ji,jj,jl) = zs2new * rswitch 456 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 457 END DO 458 END DO 459 460 ! Calculate fluxes and moments between boxes j<-->j+1 461 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 462 DO ji = fs_2, fs_jpim1 463 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 464 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt * e1v(ji,jj) / psm(ji,jj,jl) 465 zalfq = zalf * zalf 466 zalf1 = 1.0 - zalf 467 zalf1q = zalf1 * zalf1 468 ! 469 zfm (ji,jj) = zalf * psm(ji,jj,jl) 470 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 471 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 472 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 473 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 474 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 475 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 476 ! 477 ! Readjust moments remaining in the box. 478 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 479 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 480 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 481 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 482 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 483 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 484 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 485 END DO 486 END DO 487 ! 488 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 489 DO ji = fs_2, fs_jpim1 490 zalf = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) ) / psm(ji,jj+1,jl) 491 zalg (ji,jj) = zalf 492 zalfq = zalf * zalf 493 zalf1 = 1.0 - zalf 494 zalg1 (ji,jj) = zalf1 495 zalf1q = zalf1 * zalf1 496 zalg1q(ji,jj) = zalf1q 497 ! 498 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl) 499 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) & 500 & - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 501 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 502 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq 503 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 504 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl) 505 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl) 506 END DO 507 END DO 508 509 ! Readjust moments remaining in the box. 510 DO jj = 2, jpjm1 511 DO ji = fs_2, fs_jpim1 512 zbt = zbet(ji,jj-1) 513 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 514 ! 515 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 516 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 517 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 518 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 519 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 520 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 521 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 522 END DO 523 END DO 524 525 ! Put the temporary moments into appropriate neighboring boxes. 526 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 527 DO ji = fs_2, fs_jpim1 528 zbt = zbet(ji,jj-1) 529 zbt1 = 1.0 - zbet(ji,jj-1) 530 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 531 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 532 zalf1 = 1.0 - zalf 533 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 534 ! 535 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 536 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 537 & + zbt1 * psy(ji,jj,jl) 538 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 539 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 540 & + zbt1 * psyy(ji,jj,jl) 541 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 542 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 543 & + zbt1 * psxy(ji,jj,jl) 544 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 545 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 546 END DO 547 END DO 548 549 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 550 DO ji = fs_2, fs_jpim1 551 zbt = zbet(ji,jj) 552 zbt1 = 1.0 - zbet(ji,jj) 553 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 554 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 555 zalf1 = 1.0 - zalf 556 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 557 ! 558 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 559 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 560 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 561 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 562 & + ( zalf1 - zalf ) * ztemp ) ) 563 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 564 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 565 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 566 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 567 END DO 568 END DO 569 487 570 END DO 488 571 489 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)490 psm(:,:) = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )491 492 ! Calculate fluxes and moments between boxes j<-->j+1493 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0494 DO ji = 1, jpi495 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )496 zalf = MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)497 zalfq = zalf * zalf498 zalf1 = 1.0 - zalf499 zalf1q = zalf1 * zalf1500 !501 zfm (ji,jj) = zalf * psm(ji,jj)502 zf0 (ji,jj) = zalf * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj) + (zalf1-zalf) * psyy(ji,jj) ) )503 zfy (ji,jj) = zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )504 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj)505 zfx (ji,jj) = zalf * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )506 zfxy(ji,jj) = zalfq * psxy(ji,jj)507 zfxx(ji,jj) = zalf * psxx(ji,jj)508 !509 ! Readjust moments remaining in the box.510 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj)511 ps0 (ji,jj) = ps0 (ji,jj) - zf0(ji,jj)512 psy (ji,jj) = zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )513 psyy(ji,jj) = zalf1 * zalf1q * psyy(ji,jj)514 psx (ji,jj) = psx (ji,jj) - zfx(ji,jj)515 psxx(ji,jj) = psxx(ji,jj) - zfxx(ji,jj)516 psxy(ji,jj) = zalf1q * psxy(ji,jj)517 END DO518 END DO519 !520 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0.521 DO ji = 1, jpi522 zalf = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1)523 zalg (ji,jj) = zalf524 zalfq = zalf * zalf525 zalf1 = 1.0 - zalf526 zalg1 (ji,jj) = zalf1527 zalf1q = zalf1 * zalf1528 zalg1q(ji,jj) = zalf1q529 !530 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1)531 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )532 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )533 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1) * zalfq534 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )535 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1)536 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1)537 END DO538 END DO539 540 ! Readjust moments remaining in the box.541 DO jj = 2, jpj542 DO ji = 1, jpi543 zbt = zbet(ji,jj-1)544 zbt1 = ( 1.0 - zbet(ji,jj-1) )545 !546 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )547 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )548 psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )549 psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)550 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )551 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )552 psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)553 END DO554 END DO555 556 ! Put the temporary moments into appropriate neighboring boxes.557 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0.558 DO ji = 1, jpi559 zbt = zbet(ji,jj-1)560 zbt1 = ( 1.0 - zbet(ji,jj-1) )561 psm(ji,jj) = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj)562 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj)563 zalf1 = 1.0 - zalf564 ztemp = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)565 !566 ps0(ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)567 psy(ji,jj) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) &568 & + zbt1 * psy(ji,jj)569 psyy(ji,jj) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj) &570 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &571 & + zbt1 * psyy(ji,jj)572 psxy(ji,jj) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj) &573 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) ) &574 & + zbt1 * psxy(ji,jj)575 psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)576 psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)577 END DO578 END DO579 580 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0.581 DO ji = 1, jpi582 zbt = zbet(ji,jj)583 zbt1 = ( 1.0 - zbet(ji,jj) )584 psm(ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )585 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj)586 zalf1 = 1.0 - zalf587 ztemp = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)588 ps0 (ji,jj) = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )589 psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )590 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj) &591 & + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) &592 & + ( zalf1 - zalf ) * ztemp ) )593 psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) &594 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) )595 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )596 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )597 END DO598 END DO599 600 572 !-- Lateral boundary conditions 601 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T', 1., ps0 , 'T', 1. & 602 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 603 & , psxx, 'T', 1., psyy, 'T', 1. & 604 & , psxy, 'T', 1. ) 605 606 IF(ln_ctl) THEN 607 CALL prt_ctl(tab2d_1=psm , clinfo1=' adv_y: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ') 608 CALL prt_ctl(tab2d_1=psx , clinfo1=' adv_y: psx :', tab2d_2=psxx, clinfo2=' psxx : ') 609 CALL prt_ctl(tab2d_1=psy , clinfo1=' adv_y: psy :', tab2d_2=psyy, clinfo2=' psyy : ') 610 CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :') 611 ENDIF 573 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 574 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 575 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 612 576 ! 613 577 END SUBROUTINE adv_y … … 624 588 ! 625 589 ! !* allocate prather fields 626 ALLOCATE( sxopw(jpi,jpj ) , syopw(jpi,jpj) , sxxopw(jpi,jpj) , syyopw(jpi,jpj) , sxyopw(jpi,jpj), &590 ALLOCATE( sxopw(jpi,jpj,1) , syopw(jpi,jpj,1) , sxxopw(jpi,jpj,1) , syyopw(jpi,jpj,1) , sxyopw(jpi,jpj,1) , & 627 591 & sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 628 592 & sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , & … … 652 616 !! *** ROUTINE adv_pra_rst *** 653 617 !! 654 !! ** Purpose : Read or write RHGfile in restart file618 !! ** Purpose : Read or write file in restart file 655 619 !! 656 620 !! ** Method : use of IOM library -
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r10945 r11612 319 319 ! 320 320 !== Ice age ==! 321 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 322 zamsk = 1._wp 323 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 324 & poa_i, poa_i ) 325 ENDIF 321 zamsk = 1._wp 322 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 323 & poa_i, poa_i ) 326 324 ! 327 325 !== melt ponds ==!
Note: See TracChangeset
for help on using the changeset viewer.
