Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_adv_pra.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_adv_pra.F90
r10425 r11822 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 ice42 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxc0 , syc0 , sxxc0 , syyc0 , sxyc0 ! snow layers heat content 43 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxe , sye , sxxe , syye , sxye ! ice layers heat content … … 81 80 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 82 81 ! 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 82 INTEGER :: ji,jj, jk, jl, jt ! dummy loop indices 83 INTEGER :: icycle ! number of sub-timestep for the advection 84 REAL(wp) :: zdt ! - - 85 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 86 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 87 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx 88 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 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 ! transported fields 112 !------------------------- 113 z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:) ! Open water area 114 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 120 DO jk = 1, nlay_s 121 z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 122 END DO 123 DO jk = 1, nlay_i 124 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 125 END DO 126 IF ( ln_pnd_H12 ) THEN 127 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 128 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 111 ! --- transport --- ! 112 zudy(:,:) = pu_ice(:,:) * e2u(:,:) 113 zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 114 115 DO jt = 1, icycle 116 117 ! record at_i before advection (for open water) 118 zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 119 120 ! --- transported fields --- ! 121 DO jl = 1, jpl 122 zarea(:,:,jl) = e1e2t(:,:) 123 z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume 124 z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume 125 z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area 126 z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content 127 z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content 128 DO jk = 1, nlay_s 129 z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 130 END DO 131 DO jk = 1, nlay_i 132 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 133 END DO 134 IF ( ln_pnd_H12 ) THEN 135 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 136 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 137 ENDIF 138 END DO 139 ! 140 ! !--------------------------------------------! 141 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 142 ! !--------------------------------------------! 143 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 144 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 145 CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 146 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 147 CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 148 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 149 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 150 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 151 CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 152 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 153 ! 154 DO jk = 1, nlay_s !--- snow heat content 155 CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 156 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 157 CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 158 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 159 END DO 160 DO jk = 1, nlay_i !--- ice heat content 161 CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 162 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 163 CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 164 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 165 END DO 166 ! 167 IF ( ln_pnd_H12 ) THEN 168 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 169 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 170 CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 171 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 172 ENDIF 173 ! !--------------------------------------------! 174 ELSE !== even ice time step: adv_y then adv_x ==! 175 ! !--------------------------------------------! 176 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 177 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 178 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 179 CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 180 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 181 CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 182 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 183 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 184 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 185 CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 186 DO jk = 1, nlay_s !--- snow heat content 187 CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 188 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 189 CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 190 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 191 END DO 192 DO jk = 1, nlay_i !--- ice heat content 193 CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 194 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 195 CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 196 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 197 END DO 198 IF ( ln_pnd_H12 ) THEN 199 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 200 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 201 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 202 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 203 ENDIF 204 ! 129 205 ENDIF 206 207 ! --- Recover the properties from their contents --- ! 208 DO jl = 1, jpl 209 pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 210 pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 211 psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 212 poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 213 pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 214 DO jk = 1, nlay_s 215 pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 216 END DO 217 DO jk = 1, nlay_i 218 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 219 END DO 220 IF ( ln_pnd_H12 ) THEN 221 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 222 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 223 ENDIF 224 END DO 225 ! 226 ! derive open water from ice concentration 227 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 230 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water 231 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 232 END DO 233 END DO 234 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1. ) 235 ! 236 ! --- Ensure non-negative fields --- ! 237 ! Remove negative values (conservation is ensured) 238 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 239 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 ) 240 ! 241 ! --- Ensure snow load is not too big --- ! 242 CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 243 ! 130 244 END DO 131 132 ! !--------------------------------------------!133 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!134 ! !--------------------------------------------!135 DO jt = 1, initad136 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area137 & 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, jpl141 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 DO167 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 DO173 IF ( ln_pnd_H12 ) THEN174 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 ENDIF183 END DO184 END DO185 ! !--------------------------------------------!186 ELSE !== even ice time step: adv_y then adv_x ==!187 ! !--------------------------------------------!188 DO jt = 1, initad189 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area190 & 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, jpl194 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 DO220 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 DO226 IF ( ln_pnd_H12 ) THEN227 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 ENDIF236 END DO237 END DO238 ENDIF239 240 !-------------------------------------------241 ! Recover the properties from their contents242 !-------------------------------------------243 pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1)244 DO jl = 1, jpl245 pv_i (:,:, jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)246 pv_s (:,:, jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)247 psv_i(:,:, jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)248 poa_i(:,:, jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)249 pa_i (:,:, jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)250 DO jk = 1, nlay_s251 pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)252 END DO253 DO jk = 1, nlay_i254 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)255 END DO256 IF ( ln_pnd_H12 ) THEN257 pa_ip (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)258 pv_ip (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)259 ENDIF260 END DO261 !262 DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei )263 245 ! 264 246 IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file … … 267 249 268 250 269 SUBROUTINE adv_x( pd f, put , pcrh, psm , ps0 , &251 SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , & 270 252 & psx, psxx, psy , psyy, psxy ) 271 253 !!---------------------------------------------------------------------- … … 275 257 !! variable on x axis 276 258 !!---------------------------------------------------------------------- 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 moments259 REAL(wp) , INTENT(in ) :: pdt ! the time step 260 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 261 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 262 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 263 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 264 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 265 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 284 266 !! 285 INTEGER :: ji, jj 286 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! local scalars267 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 268 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 287 269 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 288 270 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 291 273 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 292 274 !----------------------------------------------------------------------- 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 275 ! 276 jcat = SIZE( ps0 , 3 ) ! size of input arrays 277 ! 278 DO jl = 1, jcat ! loop on categories 279 ! 280 ! Limitation of moments. 281 DO jj = 2, jpjm1 282 DO ji = 1, jpi 283 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 284 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 285 ! 286 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 287 zs1max = 1.5 * zslpmax 288 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 289 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 290 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 291 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 292 293 ps0 (ji,jj,jl) = zslpmax 294 psx (ji,jj,jl) = zs1new * rswitch 295 psxx(ji,jj,jl) = zs2new * rswitch 296 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 297 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 298 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 299 END DO 300 END DO 301 302 ! Calculate fluxes and moments between boxes i<-->i+1 303 DO jj = 2, jpjm1 ! Flux from i to i+1 WHEN u GT 0 304 DO ji = 1, jpi 305 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 306 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 307 zalfq = zalf * zalf 308 zalf1 = 1.0 - zalf 309 zalf1q = zalf1 * zalf1 310 ! 311 zfm (ji,jj) = zalf * psm (ji,jj,jl) 312 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 313 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 314 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 315 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 316 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 317 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 318 319 ! Readjust moments remaining in the box. 320 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 321 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 322 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 323 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 324 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 325 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 326 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 327 END DO 328 END DO 329 330 DO jj = 2, jpjm1 ! Flux from i+1 to i when u LT 0. 331 DO ji = 1, fs_jpim1 332 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 333 zalg (ji,jj) = zalf 334 zalfq = zalf * zalf 335 zalf1 = 1.0 - zalf 336 zalg1 (ji,jj) = zalf1 337 zalf1q = zalf1 * zalf1 338 zalg1q(ji,jj) = zalf1q 339 ! 340 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 341 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 342 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 343 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 344 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 345 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 346 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 347 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 348 END DO 349 END DO 350 351 DO jj = 2, jpjm1 ! Readjust moments remaining in the box. 352 DO ji = fs_2, fs_jpim1 353 zbt = zbet(ji-1,jj) 354 zbt1 = 1.0 - zbet(ji-1,jj) 355 ! 356 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 357 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 358 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 359 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 360 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 361 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 362 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 363 END DO 364 END DO 365 366 ! Put the temporary moments into appropriate neighboring boxes. 367 DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0. 368 DO ji = fs_2, fs_jpim1 369 zbt = zbet(ji-1,jj) 370 zbt1 = 1.0 - zbet(ji-1,jj) 371 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 372 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 373 zalf1 = 1.0 - zalf 374 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 375 ! 376 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 377 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 378 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 379 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 380 & + zbt1 * psxx(ji,jj,jl) 381 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 382 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 383 & + zbt1 * psxy(ji,jj,jl) 384 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 385 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 386 END DO 387 END DO 388 389 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0. 390 DO ji = fs_2, fs_jpim1 391 zbt = zbet(ji,jj) 392 zbt1 = 1.0 - zbet(ji,jj) 393 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 394 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 395 zalf1 = 1.0 - zalf 396 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 397 ! 398 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 399 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 400 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 401 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 402 & + ( zalf1 - zalf ) * ztemp ) ) 403 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 404 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 405 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 406 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 407 END DO 408 END DO 409 314 410 END DO 315 411 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 412 !-- 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 413 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 414 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 415 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 438 416 ! 439 417 END SUBROUTINE adv_x 440 418 441 419 442 SUBROUTINE adv_y( pd f, pvt , pcrh, psm , ps0 , &420 SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , & 443 421 & psx, psxx, psy , psyy, psxy ) 444 422 !!--------------------------------------------------------------------- … … 448 426 !! variable on y axis 449 427 !!--------------------------------------------------------------------- 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 moments428 REAL(wp) , INTENT(in ) :: pdt ! time step 429 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 430 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 431 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 432 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 433 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 434 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 457 435 !! 458 INTEGER :: ji, jj 459 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! temporary scalars436 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 437 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 460 438 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 461 439 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 464 442 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 465 443 !--------------------------------------------------------------------- 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 444 ! 445 jcat = SIZE( ps0 , 3 ) ! size of input arrays 446 ! 447 DO jl = 1, jcat ! loop on categories 448 ! 449 ! Limitation of moments. 450 DO jj = 1, jpj 451 DO ji = fs_2, fs_jpim1 452 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 453 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 454 ! 455 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 456 zs1max = 1.5 * zslpmax 457 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 458 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 459 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 460 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 461 ! 462 ps0 (ji,jj,jl) = zslpmax 463 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 464 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 465 psy (ji,jj,jl) = zs1new * rswitch 466 psyy(ji,jj,jl) = zs2new * rswitch 467 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 468 END DO 469 END DO 470 471 ! Calculate fluxes and moments between boxes j<-->j+1 472 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 473 DO ji = fs_2, fs_jpim1 474 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 475 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 476 zalfq = zalf * zalf 477 zalf1 = 1.0 - zalf 478 zalf1q = zalf1 * zalf1 479 ! 480 zfm (ji,jj) = zalf * psm(ji,jj,jl) 481 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 482 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 483 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 484 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 485 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 486 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 487 ! 488 ! Readjust moments remaining in the box. 489 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 490 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 491 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 492 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 493 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 494 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 495 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 496 END DO 497 END DO 498 ! 499 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 500 DO ji = fs_2, fs_jpim1 501 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 502 zalg (ji,jj) = zalf 503 zalfq = zalf * zalf 504 zalf1 = 1.0 - zalf 505 zalg1 (ji,jj) = zalf1 506 zalf1q = zalf1 * zalf1 507 zalg1q(ji,jj) = zalf1q 508 ! 509 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl) 510 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) & 511 & - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 512 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 513 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq 514 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 515 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl) 516 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl) 517 END DO 518 END DO 519 520 ! Readjust moments remaining in the box. 521 DO jj = 2, jpjm1 522 DO ji = fs_2, fs_jpim1 523 zbt = zbet(ji,jj-1) 524 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 525 ! 526 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 527 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 528 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 529 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 530 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 531 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 532 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 533 END DO 534 END DO 535 536 ! Put the temporary moments into appropriate neighboring boxes. 537 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 538 DO ji = fs_2, fs_jpim1 539 zbt = zbet(ji,jj-1) 540 zbt1 = 1.0 - zbet(ji,jj-1) 541 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 542 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 543 zalf1 = 1.0 - zalf 544 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 545 ! 546 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 547 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 548 & + zbt1 * psy(ji,jj,jl) 549 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 550 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 551 & + zbt1 * psyy(ji,jj,jl) 552 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 553 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 554 & + zbt1 * psxy(ji,jj,jl) 555 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 556 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 557 END DO 558 END DO 559 560 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 561 DO ji = fs_2, fs_jpim1 562 zbt = zbet(ji,jj) 563 zbt1 = 1.0 - zbet(ji,jj) 564 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 565 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 566 zalf1 = 1.0 - zalf 567 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 568 ! 569 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 570 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 571 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 572 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 573 & + ( zalf1 - zalf ) * ztemp ) ) 574 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 575 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 576 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 577 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 578 END DO 579 END DO 580 487 581 END DO 488 582 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) 583 !-- Lateral boundary conditions 584 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 585 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 586 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 587 ! 588 END SUBROUTINE adv_y 589 590 591 SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 592 !!------------------------------------------------------------------- 593 !! *** ROUTINE Hsnow *** 594 !! 595 !! ** Purpose : 1- Check snow load after advection 596 !! 2- Correct pond concentration to avoid a_ip > a_i 597 !! 598 !! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface 599 !! then put the snow excess in the ocean 600 !! 601 !! ** Notes : This correction is crucial because of the call to routine icecor afterwards 602 !! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 603 !! make the snow very thick (if concentration decreases drastically) 604 !! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 605 !!------------------------------------------------------------------- 606 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 607 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip 608 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 609 ! 610 INTEGER :: ji, jj, jl ! dummy loop indices 611 REAL(wp) :: z1_dt, zvs_excess, zfra 612 !!------------------------------------------------------------------- 613 ! 614 z1_dt = 1._wp / pdt 615 ! 616 ! -- check snow load -- ! 617 DO jl = 1, jpl 618 DO jj = 1, jpj 619 DO ji = 1, jpi 620 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 621 ! 622 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 623 ! 624 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 625 ! put snow excess in the ocean 626 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 627 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 628 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 629 ! correct snow volume and heat content 630 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 631 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 632 ENDIF 633 ! 634 ENDIF 635 END DO 517 636 END DO 518 637 END DO 519 638 ! 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 639 !-- correct pond concentration to avoid a_ip > a_i -- ! 640 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 641 ! 642 END SUBROUTINE Hsnow 614 643 615 644 … … 624 653 ! 625 654 ! !* allocate prather fields 626 ALLOCATE( sxopw(jpi,jpj) , syopw(jpi,jpj) , sxxopw(jpi,jpj) , syyopw(jpi,jpj) , sxyopw(jpi,jpj) , & 627 & sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 655 ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 628 656 & sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , & 629 657 & sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , & … … 652 680 !! *** ROUTINE adv_pra_rst *** 653 681 !! 654 !! ** Purpose : Read or write RHGfile in restart file682 !! ** Purpose : Read or write file in restart file 655 683 !! 656 684 !! ** Method : use of IOM library … … 671 699 ! !==========================! 672 700 ! 673 IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sx opw' , ldstop = .FALSE. ) ! file exist: id1>0701 IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. ) ! file exist: id1>0 674 702 ELSE ; id1 = 0 ! no restart: id1=0 675 703 ENDIF … … 689 717 CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn ) 690 718 CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn ) 691 ! ! lead fraction719 ! ! ice concentration 692 720 CALL iom_get( numrir, jpdom_autoglo, 'sxa' , sxa ) 693 721 CALL iom_get( numrir, jpdom_autoglo, 'sya' , sya ) … … 707 735 CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 708 736 CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage ) 709 ! ! open water in sea ice710 CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw )711 CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw )712 CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )713 CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )714 CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )715 737 ! ! snow layers heat content 716 738 DO jk = 1, nlay_s … … 752 774 sxice = 0._wp ; syice = 0._wp ; sxxice = 0._wp ; syyice = 0._wp ; sxyice = 0._wp ! ice thickness 753 775 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 fraction776 sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! ice concentration 755 777 sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity 756 778 sxage = 0._wp ; syage = 0._wp ; sxxage = 0._wp ; syyage = 0._wp ; sxyage = 0._wp ! ice age 757 sxopw = 0._wp ; syopw = 0._wp ; sxxopw = 0._wp ; syyopw = 0._wp ; sxyopw = 0._wp ! open water in sea ice758 779 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 759 780 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content … … 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 ) … … 804 825 CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 805 826 CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage ) 806 ! ! open water in sea ice807 CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw )808 CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw )809 CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )810 CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )811 CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )812 827 ! ! snow layers heat content 813 828 DO jk = 1, nlay_s
Note: See TracChangeset
for help on using the changeset viewer.