Changeset 11732
- Timestamp:
- 2019-10-18T18:10:30+02:00 (5 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 4 edited
-
icedyn_adv_pra.F90 (modified) (12 diffs)
-
icedyn_rdgrft.F90 (modified) (11 diffs)
-
iceitd.F90 (modified) (4 diffs)
-
icevar.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv_pra.F90
r11632 r11732 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 :: j k, jl, jt! dummy loop indices82 INTEGER :: ji,jj, jk, jl, jt ! dummy loop indices 84 83 INTEGER :: icycle ! number of sub-timestep for the advection 85 84 REAL(wp) :: zdt ! - - 86 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 87 88 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 88 REAL(wp), DIMENSION(jpi,jpj,1) :: z0opw89 89 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 90 90 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp … … 109 109 zdt = rdt_ice / REAL(icycle) 110 110 111 !------------------------- 112 ! transported fields 113 !------------------------- 114 z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:) ! Open water area 115 DO jl = 1, jpl 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 122 DO jk = 1, nlay_s 123 z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 124 END DO 125 DO jk = 1, nlay_i 126 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 127 END DO 128 IF ( ln_pnd_H12 ) THEN 129 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 130 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 131 ENDIF 132 END DO 133 134 ! !--------------------------------------------! 135 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 136 ! !--------------------------------------------! 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 ) 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 ) 150 153 ! 151 DO jk = 1, nlay_s 152 CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &153 & 154 CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &155 & 156 END DO 157 DO jk = 1, nlay_i 158 CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &159 & 160 CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &161 & 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,:) ) 162 165 END DO 163 166 ! 164 167 IF ( ln_pnd_H12 ) THEN 165 CALL adv_x( zdt , pu_ice, 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction166 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 volume168 CALL adv_y( zdt , pv_ice, 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )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 ) 169 172 ENDIF 170 END DO 171 ! !--------------------------------------------! 172 ELSE !== even ice time step: adv_y then adv_x ==! 173 ! !--------------------------------------------! 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,:) ) 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,:) ) 198 197 END DO 199 198 IF ( ln_pnd_H12 ) THEN 200 CALL adv_y( zdt , pv_ice, 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction201 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 volume203 CALL adv_x( zdt , pu_ice, 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )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 ) 204 203 ENDIF 205 END DO 206 ENDIF 207 208 !------------------------------------------- 209 ! Recover the properties from their contents 210 !------------------------------------------- 211 pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1) 212 DO jl = 1, jpl 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) 218 DO jk = 1, nlay_s 219 pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 220 END DO 221 DO jk = 1, nlay_i 222 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 223 END DO 224 IF ( ln_pnd_H12 ) THEN 225 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 226 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 204 ! 227 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 ! 228 244 END DO 229 !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 )237 245 ! 238 246 IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file … … 296 304 DO ji = 1, jpi 297 305 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)306 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 299 307 zalfq = zalf * zalf 300 308 zalf1 = 1.0 - zalf … … 322 330 DO jj = 2, jpjm1 ! Flux from i+1 to i when u LT 0. 323 331 DO ji = 1, fs_jpim1 324 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj)/ psm(ji+1,jj,jl)332 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 325 333 zalg (ji,jj) = zalf 326 334 zalfq = zalf * zalf … … 465 473 DO ji = fs_2, fs_jpim1 466 474 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)475 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 468 476 zalfq = zalf * zalf 469 477 zalf1 = 1.0 - zalf … … 491 499 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 492 500 DO ji = fs_2, fs_jpim1 493 zalf = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) )/ psm(ji,jj+1,jl)501 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 494 502 zalg (ji,jj) = zalf 495 503 zalfq = zalf * zalf … … 645 653 ! 646 654 ! !* allocate prather fields 647 ALLOCATE( sxopw(jpi,jpj,1) , syopw(jpi,jpj,1) , sxxopw(jpi,jpj,1) , syyopw(jpi,jpj,1) , sxyopw(jpi,jpj,1) , & 648 & 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) , & 649 656 & sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , & 650 657 & sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , & … … 692 699 ! !==========================! 693 700 ! 694 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 695 702 ELSE ; id1 = 0 ! no restart: id1=0 696 703 ENDIF … … 728 735 CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 729 736 CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage ) 730 ! ! open water in sea ice731 CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw )732 CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw )733 CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )734 CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )735 CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )736 737 ! ! snow layers heat content 737 738 DO jk = 1, nlay_s … … 776 777 sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity 777 778 sxage = 0._wp ; syage = 0._wp ; sxxage = 0._wp ; syyage = 0._wp ; sxyage = 0._wp ! ice age 778 sxopw = 0._wp ; syopw = 0._wp ; sxxopw = 0._wp ; syyopw = 0._wp ; sxyopw = 0._wp ! open water in sea ice779 779 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 780 780 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content … … 825 825 CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 826 826 CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage ) 827 ! ! open water in sea ice828 CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw )829 CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw )830 CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )831 CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )832 CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )833 827 ! ! snow layers heat content 834 828 DO jk = 1, nlay_s -
NEMO/trunk/src/ICE/icedyn_rdgrft.F90
r11560 r11732 86 86 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 87 !!------------------------------------------------------------------- 88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij),&89 & apartf(jpij,0:jpl) , hrmin(jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl),&90 & hrmax (jpij,jpl), hi_hrdg(jpij,jpl) , araft (jpij,jpl),&88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 89 & apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 90 & hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , & 91 91 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 92 92 … … 137 137 REAL(wp) :: zfac ! local scalar 138 138 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 139 REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s)140 139 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 141 140 ! … … 175 174 176 175 ! just needed here 177 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i )178 176 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 179 177 ! needed here and in the iteration loop 178 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_i) 180 179 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 181 180 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) … … 187 186 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 188 187 ! 189 ! divergence given by the advection scheme 190 ! (which may not be equal to divu as computed from the velocity field) 191 IF ( ln_adv_Pra ) THEN 192 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 193 ELSEIF( ln_adv_UMx ) THEN 194 zdivu_adv(ji) = zdivu(ji) 195 ENDIF 196 ! 197 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! make sure the closing rate is large enough 198 ! ! to give asum = 1.0 after ridging 188 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough 189 ! ! to give asum = 1.0 after ridging 199 190 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 200 opning(ji) = closing_net(ji) + zdivu _adv(ji)191 opning(ji) = closing_net(ji) + zdivu(ji) 201 192 END DO 202 193 ! … … 215 206 ato_i_1d (ipti) = ato_i_1d (ji) 216 207 closing_net(ipti) = closing_net(ji) 217 zdivu _adv (ipti) = zdivu_adv(ji)208 zdivu (ipti) = zdivu (ji) 218 209 opning (ipti) = opning (ji) 219 210 ENDIF … … 259 250 ELSE 260 251 iterate_ridging = 1 261 zdivu _adv(ji) = zfac * r1_rdtice262 closing_net(ji) = MAX( 0._wp, -zdivu _adv(ji) )263 opning (ji) = MAX( 0._wp, zdivu _adv(ji) )252 zdivu (ji) = zfac * r1_rdtice 253 closing_net(ji) = MAX( 0._wp, -zdivu(ji) ) 254 opning (ji) = MAX( 0._wp, zdivu(ji) ) 264 255 ENDIF 265 256 END DO … … 309 300 310 301 ! ! Ice thickness needed for rafting 311 WHERE( pa_i(1:npti,:) > epsi 20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)302 WHERE( pa_i(1:npti,:) > epsi10 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 312 303 ELSEWHERE ; zhi(1:npti,:) = 0._wp 313 304 END WHERE … … 328 319 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 329 320 ! 330 WHERE( zasum(1:npti) > epsi 20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)321 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 331 322 ELSEWHERE ; z1_asum(1:npti) = 0._wp 332 323 END WHERE … … 454 445 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 455 446 ! NOTE: 0 < aksum <= 1 456 WHERE( zaksum(1:npti) > epsi 20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)447 WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 457 448 ELSEWHERE ; closing_gross(1:npti) = 0._wp 458 449 END WHERE … … 537 528 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 538 529 539 IF( a_i_2d(ji,jl1) > epsi 20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)530 IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 540 531 ELSE ; z1_ai(ji) = 0._wp 541 532 ENDIF … … 595 586 ! virtual salt flux to keep salinity constant 596 587 IF( nn_icesal /= 2 ) THEN 597 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) 588 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i 598 589 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice & ! put back sss_m into the ocean 599 590 & - s_i_1d(ji) * vsw * rhoi * r1_rdtice ! and get s_i from the ocean -
NEMO/trunk/src/ICE/iceitd.F90
r11536 r11732 211 211 CALL itd_glinear( zhb0(1:npti) , zhb1(1:npti) , h_ib_1d(1:npti) , a_i_1d(1:npti) , & ! in 212 212 & g0 (1:npti,1), g1 (1:npti,1), hL (1:npti,1), hR (1:npti,1) ) ! out 213 213 ! 214 214 ! Area lost due to melting of thin ice 215 215 DO ji = 1, npti … … 218 218 ! 219 219 zdh0 = h_i_1d(ji) - h_ib_1d(ji) 220 IF( zdh0 < 0.0 ) THEN ! remove area from category 1220 IF( zdh0 < 0.0 ) THEN ! remove area from category 1 221 221 zdh0 = MIN( -zdh0, hi_max(1) ) 222 222 !Integrate g(1) from 0 to dh0 to estimate area melted … … 226 226 zx1 = zetamax 227 227 zx2 = 0.5 * zetamax * zetamax 228 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 228 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed 229 229 zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i 230 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 231 ! of thin ice (zdamax > 0) 230 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting of thin ice (zdamax > 0) 232 231 ! Remove area, conserving volume 233 232 h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) … … 349 348 DO ji = 1, npti 350 349 ! 351 IF( paice(ji) > epsi10 .AND. phice(ji) > 0._wp) THEN350 IF( paice(ji) > epsi10 .AND. phice(ji) > epsi10 ) THEN 352 351 ! 353 352 ! Initialize hL and hR -
NEMO/trunk/src/ICE/icevar.F90
r11560 r11732 622 622 pv_s (ji,jj,jl) = 0._wp 623 623 ENDIF 624 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN624 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN 625 625 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 626 626 psv_i (ji,jj,jl) = 0._wp
Note: See TracChangeset
for help on using the changeset viewer.
