Changeset 11692 for NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/ICE/icedyn_adv_pra.F90
- Timestamp:
- 2019-10-12T16:08:18+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/ICE/icedyn_adv_pra.F90
r10425 r11692 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 … … 36 36 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxice, syice, sxxice, syyice, sxyice ! ice thickness 37 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsn , sysn , sxxsn , syysn , sxysn ! snow thickness 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! lead fraction38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! ice concentration 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 ) 234 ! 235 ! --- Ensure snow load is not too big --- ! 236 CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 263 237 ! 264 238 IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file … … 267 241 268 242 269 SUBROUTINE adv_x( pd f, put , pcrh, psm , ps0 , &243 SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , & 270 244 & psx, psxx, psy , psyy, psxy ) 271 245 !!---------------------------------------------------------------------- … … 275 249 !! variable on x axis 276 250 !!---------------------------------------------------------------------- 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 moments251 REAL(wp) , INTENT(in ) :: pdt ! the time step 252 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 253 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 254 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 255 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 256 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 257 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 284 258 !! 285 INTEGER :: ji, jj 286 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! local scalars259 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 260 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 287 261 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 288 262 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 291 265 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 292 266 !----------------------------------------------------------------------- 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 267 ! 268 jcat = SIZE( ps0 , 3 ) ! size of input arrays 269 ! 270 DO jl = 1, jcat ! loop on categories 271 ! 272 ! Limitation of moments. 273 DO jj = 2, jpjm1 274 DO ji = 1, jpi 275 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 276 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 277 ! 278 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 279 zs1max = 1.5 * zslpmax 280 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 281 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 282 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 283 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 284 285 ps0 (ji,jj,jl) = zslpmax 286 psx (ji,jj,jl) = zs1new * rswitch 287 psxx(ji,jj,jl) = zs2new * rswitch 288 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 289 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 290 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 291 END DO 292 END DO 293 294 ! Calculate fluxes and moments between boxes i<-->i+1 295 DO jj = 2, jpjm1 ! Flux from i to i+1 WHEN u GT 0 296 DO ji = 1, jpi 297 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 298 zalf = MAX( 0._wp, put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji,jj,jl) 299 zalfq = zalf * zalf 300 zalf1 = 1.0 - zalf 301 zalf1q = zalf1 * zalf1 302 ! 303 zfm (ji,jj) = zalf * psm (ji,jj,jl) 304 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 305 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 306 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 307 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 308 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 309 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 310 311 ! Readjust moments remaining in the box. 312 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 313 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 314 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 315 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 316 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 317 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 318 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 319 END DO 320 END DO 321 322 DO jj = 2, jpjm1 ! Flux from i+1 to i when u LT 0. 323 DO ji = 1, fs_jpim1 324 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji+1,jj,jl) 325 zalg (ji,jj) = zalf 326 zalfq = zalf * zalf 327 zalf1 = 1.0 - zalf 328 zalg1 (ji,jj) = zalf1 329 zalf1q = zalf1 * zalf1 330 zalg1q(ji,jj) = zalf1q 331 ! 332 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 333 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 334 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 335 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 336 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 337 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 338 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 339 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 340 END DO 341 END DO 342 343 DO jj = 2, jpjm1 ! Readjust moments remaining in the box. 344 DO ji = fs_2, fs_jpim1 345 zbt = zbet(ji-1,jj) 346 zbt1 = 1.0 - zbet(ji-1,jj) 347 ! 348 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 349 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 350 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 351 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 352 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 353 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 354 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 355 END DO 356 END DO 357 358 ! Put the temporary moments into appropriate neighboring boxes. 359 DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0. 360 DO ji = fs_2, fs_jpim1 361 zbt = zbet(ji-1,jj) 362 zbt1 = 1.0 - zbet(ji-1,jj) 363 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 364 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 365 zalf1 = 1.0 - zalf 366 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 367 ! 368 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 369 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 370 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 371 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 372 & + zbt1 * psxx(ji,jj,jl) 373 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 374 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 375 & + zbt1 * psxy(ji,jj,jl) 376 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 377 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 378 END DO 379 END DO 380 381 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0. 382 DO ji = fs_2, fs_jpim1 383 zbt = zbet(ji,jj) 384 zbt1 = 1.0 - zbet(ji,jj) 385 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 386 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 387 zalf1 = 1.0 - zalf 388 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 389 ! 390 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 391 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 392 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 393 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 394 & + ( zalf1 - zalf ) * ztemp ) ) 395 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 396 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 397 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 398 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 399 END DO 400 END DO 401 314 402 END DO 315 403 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 404 !-- 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 405 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 406 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 407 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 438 408 ! 439 409 END SUBROUTINE adv_x 440 410 441 411 442 SUBROUTINE adv_y( pd f, pvt , pcrh, psm , ps0 , &412 SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , & 443 413 & psx, psxx, psy , psyy, psxy ) 444 414 !!--------------------------------------------------------------------- … … 448 418 !! variable on y axis 449 419 !!--------------------------------------------------------------------- 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 moments420 REAL(wp) , INTENT(in ) :: pdt ! time step 421 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 422 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 423 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 424 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 425 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 426 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 457 427 !! 458 INTEGER :: ji, jj 459 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! temporary scalars428 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 429 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 460 430 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 461 431 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 464 434 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 465 435 !--------------------------------------------------------------------- 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 436 ! 437 jcat = SIZE( ps0 , 3 ) ! size of input arrays 438 ! 439 DO jl = 1, jcat ! loop on categories 440 ! 441 ! Limitation of moments. 442 DO jj = 1, jpj 443 DO ji = fs_2, fs_jpim1 444 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 445 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 446 ! 447 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 448 zs1max = 1.5 * zslpmax 449 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 450 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 451 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 452 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 453 ! 454 ps0 (ji,jj,jl) = zslpmax 455 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 456 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 457 psy (ji,jj,jl) = zs1new * rswitch 458 psyy(ji,jj,jl) = zs2new * rswitch 459 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 460 END DO 461 END DO 462 463 ! Calculate fluxes and moments between boxes j<-->j+1 464 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 465 DO ji = fs_2, fs_jpim1 466 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 467 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt * e1v(ji,jj) / psm(ji,jj,jl) 468 zalfq = zalf * zalf 469 zalf1 = 1.0 - zalf 470 zalf1q = zalf1 * zalf1 471 ! 472 zfm (ji,jj) = zalf * psm(ji,jj,jl) 473 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 474 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 475 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 476 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 477 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 478 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 479 ! 480 ! Readjust moments remaining in the box. 481 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 482 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 483 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 484 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 485 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 486 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 487 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 488 END DO 489 END DO 490 ! 491 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 492 DO ji = fs_2, fs_jpim1 493 zalf = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) ) / psm(ji,jj+1,jl) 494 zalg (ji,jj) = zalf 495 zalfq = zalf * zalf 496 zalf1 = 1.0 - zalf 497 zalg1 (ji,jj) = zalf1 498 zalf1q = zalf1 * zalf1 499 zalg1q(ji,jj) = zalf1q 500 ! 501 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl) 502 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) & 503 & - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 504 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 505 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq 506 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 507 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl) 508 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl) 509 END DO 510 END DO 511 512 ! Readjust moments remaining in the box. 513 DO jj = 2, jpjm1 514 DO ji = fs_2, fs_jpim1 515 zbt = zbet(ji,jj-1) 516 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 517 ! 518 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 519 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 520 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 521 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 522 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 523 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 524 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 525 END DO 526 END DO 527 528 ! Put the temporary moments into appropriate neighboring boxes. 529 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 530 DO ji = fs_2, fs_jpim1 531 zbt = zbet(ji,jj-1) 532 zbt1 = 1.0 - zbet(ji,jj-1) 533 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 534 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 535 zalf1 = 1.0 - zalf 536 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 537 ! 538 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 539 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 540 & + zbt1 * psy(ji,jj,jl) 541 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 542 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 543 & + zbt1 * psyy(ji,jj,jl) 544 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 545 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 546 & + zbt1 * psxy(ji,jj,jl) 547 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 548 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 549 END DO 550 END DO 551 552 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 553 DO ji = fs_2, fs_jpim1 554 zbt = zbet(ji,jj) 555 zbt1 = 1.0 - zbet(ji,jj) 556 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 557 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 558 zalf1 = 1.0 - zalf 559 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 560 ! 561 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 562 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 563 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 564 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 565 & + ( zalf1 - zalf ) * ztemp ) ) 566 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 567 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 568 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 569 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 570 END DO 571 END DO 572 487 573 END DO 488 574 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+1 493 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 494 DO ji = 1, jpi 495 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 * zalf 498 zalf1 = 1.0 - zalf 499 zalf1q = zalf1 * zalf1 500 ! 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) 575 !-- Lateral boundary conditions 576 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 577 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 578 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 579 ! 580 END SUBROUTINE adv_y 581 582 583 SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 584 !!------------------------------------------------------------------- 585 !! *** ROUTINE Hsnow *** 586 !! 587 !! ** Purpose : 1- Check snow load after advection 588 !! 2- Correct pond concentration to avoid a_ip > a_i 589 !! 590 !! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface 591 !! then put the snow excess in the ocean 592 !! 593 !! ** Notes : This correction is crucial because of the call to routine icecor afterwards 594 !! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 595 !! make the snow very thick (if concentration decreases drastically) 596 !! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 597 !!------------------------------------------------------------------- 598 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 599 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip 600 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 601 ! 602 INTEGER :: ji, jj, jl ! dummy loop indices 603 REAL(wp) :: z1_dt, zvs_excess, zfra 604 !!------------------------------------------------------------------- 605 ! 606 z1_dt = 1._wp / pdt 607 ! 608 ! -- check snow load -- ! 609 DO jl = 1, jpl 610 DO jj = 1, jpj 611 DO ji = 1, jpi 612 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 613 ! 614 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 615 ! 616 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 617 ! put snow excess in the ocean 618 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 619 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 620 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 621 ! correct snow volume and heat content 622 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 623 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 624 ENDIF 625 ! 626 ENDIF 627 END DO 517 628 END DO 518 629 END DO 519 630 ! 520 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 521 DO ji = 1, jpi 522 zalf = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 523 zalg (ji,jj) = zalf 524 zalfq = zalf * zalf 525 zalf1 = 1.0 - zalf 526 zalg1 (ji,jj) = zalf1 527 zalf1q = zalf1 * zalf1 528 zalg1q(ji,jj) = zalf1q 529 ! 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) * zalfq 534 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 DO 538 END DO 539 540 ! Readjust moments remaining in the box. 541 DO jj = 2, jpj 542 DO ji = 1, jpi 543 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 DO 554 END DO 555 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, jpi 559 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 - zalf 564 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 DO 578 END DO 579 580 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 581 DO ji = 1, jpi 582 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 - zalf 587 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 DO 598 END DO 599 600 !-- 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 612 ! 613 END SUBROUTINE adv_y 631 !-- correct pond concentration to avoid a_ip > a_i -- ! 632 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 633 ! 634 END SUBROUTINE Hsnow 614 635 615 636 … … 624 645 ! 625 646 ! !* allocate prather fields 626 ALLOCATE( sxopw(jpi,jpj ) , syopw(jpi,jpj) , sxxopw(jpi,jpj) , syyopw(jpi,jpj) , sxyopw(jpi,jpj), &647 ALLOCATE( sxopw(jpi,jpj,1) , syopw(jpi,jpj,1) , sxxopw(jpi,jpj,1) , syyopw(jpi,jpj,1) , sxyopw(jpi,jpj,1) , & 627 648 & sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 628 649 & sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , & … … 652 673 !! *** ROUTINE adv_pra_rst *** 653 674 !! 654 !! ** Purpose : Read or write RHGfile in restart file675 !! ** Purpose : Read or write file in restart file 655 676 !! 656 677 !! ** Method : use of IOM library … … 689 710 CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn ) 690 711 CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn ) 691 ! ! lead fraction712 ! ! ice concentration 692 713 CALL iom_get( numrir, jpdom_autoglo, 'sxa' , sxa ) 693 714 CALL iom_get( numrir, jpdom_autoglo, 'sya' , sya ) … … 752 773 sxice = 0._wp ; syice = 0._wp ; sxxice = 0._wp ; syyice = 0._wp ; sxyice = 0._wp ! ice thickness 753 774 sxsn = 0._wp ; sysn = 0._wp ; sxxsn = 0._wp ; syysn = 0._wp ; sxysn = 0._wp ! snow thickness 754 sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! lead fraction775 sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! ice concentration 755 776 sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity 756 777 sxage = 0._wp ; syage = 0._wp ; sxxage = 0._wp ; syyage = 0._wp ; sxyage = 0._wp ! ice age … … 786 807 CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn ) 787 808 CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn ) 788 ! ! lead fraction809 ! ! ice concentration 789 810 CALL iom_rstput( iter, nitrst, numriw, 'sxa' , sxa ) 790 811 CALL iom_rstput( iter, nitrst, numriw, 'sya' , sya )
Note: See TracChangeset
for help on using the changeset viewer.
