- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_pra.F90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icedyn_adv_pra.F90
r12178 r12928 16 16 !! adv_pra_rst : read/write Prather field in ice restart file, or initialized to zero 17 17 !!---------------------------------------------------------------------- 18 USE phycst ! physical constant 18 19 USE dom_oce ! ocean domain 19 20 USE ice ! sea-ice variables … … 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsal, sysal, sxxsal, syysal, sxysal ! ice salinity 40 41 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 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 … … 46 46 47 47 !! * Substitutions 48 # include " vectopt_loop_substitute.h90"48 # include "do_loop_substitute.h90" 49 49 !!---------------------------------------------------------------------- 50 50 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 54 54 CONTAINS 55 55 56 SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, &56 SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 57 57 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 58 58 !!---------------------------------------------------------------------- … … 70 70 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity 71 71 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity 72 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness 73 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness 74 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness 72 75 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 73 76 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 81 84 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 82 85 ! 83 INTEGER :: j k, jl, jt! dummy loop indices86 INTEGER :: ji,jj, jk, jl, jt ! dummy loop indices 84 87 INTEGER :: icycle ! number of sub-timestep for the advection 85 88 REAL(wp) :: zdt ! - - 86 89 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 90 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 91 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx 92 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 87 93 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 88 REAL(wp), DIMENSION(jpi,jpj,1) :: z0opw89 94 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 90 95 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp … … 95 100 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 96 101 ! 102 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 103 DO jl = 1, jpl 104 DO_2D_00_00 105 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 106 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 107 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 108 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 109 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 110 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 111 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 112 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 113 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 114 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 115 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 116 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 117 END_2D 118 END DO 119 CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 120 ! 97 121 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! 98 122 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 99 123 ! this should not affect too much the stability 100 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * r dt_ice * r1_e1u(:,:) )101 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * r dt_ice * r1_e2v(:,:) ) )124 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rDt_ice * r1_e1u(:,:) ) 125 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rDt_ice * r1_e2v(:,:) ) ) 102 126 103 127 ! non-blocking global communication send zcflnow and receive zcflprv … … 107 131 ELSE ; icycle = 1 108 132 ENDIF 109 zdt = r dt_ice / REAL(icycle)133 zdt = rDt_ice / REAL(icycle) 110 134 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 ) 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,:) ) 135 ! --- transport --- ! 136 zudy(:,:) = pu_ice(:,:) * e2u(:,:) 137 zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 138 139 DO jt = 1, icycle 140 141 ! record at_i before advection (for open water) 142 zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 143 144 ! --- transported fields --- ! 145 DO jl = 1, jpl 146 zarea(:,:,jl) = e1e2t(:,:) 147 z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume 148 z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume 149 z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area 150 z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content 151 z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content 152 DO jk = 1, nlay_s 153 z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 156 154 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,:) ) 155 DO jk = 1, nlay_i 156 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 162 157 END DO 163 !164 158 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 ) 159 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 160 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 169 161 ENDIF 170 162 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,:) ) 163 ! 164 ! !--------------------------------------------! 165 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 166 ! !--------------------------------------------! 167 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 168 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 169 CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 170 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 171 CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 172 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 173 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 174 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 175 CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 176 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 177 ! 178 DO jk = 1, nlay_s !--- snow heat content 179 CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 180 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 181 CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 182 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 192 183 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,:) ) 184 DO jk = 1, nlay_i !--- ice heat content 185 CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 186 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 187 CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 188 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 189 END DO 190 ! 191 IF ( ln_pnd_H12 ) THEN 192 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 193 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 194 CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 195 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 196 ENDIF 197 ! !--------------------------------------------! 198 ELSE !== even ice time step: adv_y then adv_x ==! 199 ! !--------------------------------------------! 200 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 201 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 202 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 203 CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 204 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 205 CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 206 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 207 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 208 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 209 CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 210 DO jk = 1, nlay_s !--- snow heat content 211 CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 212 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 213 CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 214 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 215 END DO 216 DO jk = 1, nlay_i !--- ice heat content 217 CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 218 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 219 CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 220 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 198 221 END DO 199 222 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 ) 223 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 224 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 225 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 226 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 227 ENDIF 228 ! 229 ENDIF 230 231 ! --- Recover the properties from their contents --- ! 232 DO jl = 1, jpl 233 pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 234 pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 235 psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 236 poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 237 pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 238 DO jk = 1, nlay_s 239 pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 240 END DO 241 DO jk = 1, nlay_i 242 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 243 END DO 244 IF ( ln_pnd_H12 ) THEN 245 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 246 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 204 247 ENDIF 205 248 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) 227 ENDIF 249 ! 250 ! derive open water from ice concentration 251 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 252 DO_2D_00_00 253 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water 254 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 255 END_2D 256 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1. ) 257 ! 258 ! --- Ensure non-negative fields --- ! 259 ! Remove negative values (conservation is ensured) 260 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 261 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 ) 262 ! 263 ! --- Make sure ice thickness is not too big --- ! 264 ! (because ice thickness can be too large where ice concentration is very small) 265 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 266 ! 267 ! --- Ensure snow load is not too big --- ! 268 CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 269 ! 228 270 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 271 ! 238 272 IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file … … 271 305 ! 272 306 ! 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 307 DO_2D_00_11 308 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 309 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 310 ! 311 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 312 zs1max = 1.5 * zslpmax 313 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 314 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 315 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 316 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 317 318 ps0 (ji,jj,jl) = zslpmax 319 psx (ji,jj,jl) = zs1new * rswitch 320 psxx(ji,jj,jl) = zs2new * rswitch 321 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 322 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 323 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 324 END_2D 293 325 294 326 ! 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 327 DO_2D_00_11 328 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 329 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 330 zalfq = zalf * zalf 331 zalf1 = 1.0 - zalf 332 zalf1q = zalf1 * zalf1 333 ! 334 zfm (ji,jj) = zalf * psm (ji,jj,jl) 335 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 336 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 337 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 338 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 339 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 340 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 341 342 ! Readjust moments remaining in the box. 343 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 344 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 345 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 346 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 347 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 348 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 349 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 350 END_2D 351 352 DO_2D_00_10 353 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 354 zalg (ji,jj) = zalf 355 zalfq = zalf * zalf 356 zalf1 = 1.0 - zalf 357 zalg1 (ji,jj) = zalf1 358 zalf1q = zalf1 * zalf1 359 zalg1q(ji,jj) = zalf1q 360 ! 361 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 362 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 363 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 364 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 365 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 366 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 367 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 368 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 369 END_2D 370 371 DO_2D_00_00 372 zbt = zbet(ji-1,jj) 373 zbt1 = 1.0 - zbet(ji-1,jj) 374 ! 375 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 376 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 377 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 378 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 379 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 380 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 381 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 382 END_2D 357 383 358 384 ! 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 385 DO_2D_00_00 386 zbt = zbet(ji-1,jj) 387 zbt1 = 1.0 - zbet(ji-1,jj) 388 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 389 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 390 zalf1 = 1.0 - zalf 391 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 392 ! 393 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 394 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 395 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 396 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 397 & + zbt1 * psxx(ji,jj,jl) 398 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 399 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 400 & + zbt1 * psxy(ji,jj,jl) 401 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 402 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 403 END_2D 404 405 DO_2D_00_00 406 zbt = zbet(ji,jj) 407 zbt1 = 1.0 - zbet(ji,jj) 408 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 409 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 410 zalf1 = 1.0 - zalf 411 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 412 ! 413 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 414 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 415 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 416 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 417 & + ( zalf1 - zalf ) * ztemp ) ) 418 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 419 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 420 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 421 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 422 END_2D 401 423 402 424 END DO … … 440 462 ! 441 463 ! 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 464 DO_2D_11_00 465 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 466 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 467 ! 468 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 469 zs1max = 1.5 * zslpmax 470 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 471 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 472 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 473 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 474 ! 475 ps0 (ji,jj,jl) = zslpmax 476 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 477 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 478 psy (ji,jj,jl) = zs1new * rswitch 479 psyy(ji,jj,jl) = zs2new * rswitch 480 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 481 END_2D 462 482 463 483 ! 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 484 DO_2D_11_00 485 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 486 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 487 zalfq = zalf * zalf 488 zalf1 = 1.0 - zalf 489 zalf1q = zalf1 * zalf1 490 ! 491 zfm (ji,jj) = zalf * psm(ji,jj,jl) 492 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 493 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 494 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 495 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 496 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 497 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 498 ! 499 ! Readjust moments remaining in the box. 500 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 501 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 502 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 503 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 504 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 505 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 506 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 507 END_2D 508 ! 509 DO_2D_10_00 510 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 511 zalg (ji,jj) = zalf 512 zalfq = zalf * zalf 513 zalf1 = 1.0 - zalf 514 zalg1 (ji,jj) = zalf1 515 zalf1q = zalf1 * zalf1 516 zalg1q(ji,jj) = zalf1q 517 ! 518 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl) 519 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) & 520 & - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 521 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 522 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq 523 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 524 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl) 525 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl) 526 END_2D 511 527 512 528 ! 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 529 DO_2D_00_00 530 zbt = zbet(ji,jj-1) 531 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 532 ! 533 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 534 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 535 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 536 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 537 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 538 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 539 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 540 END_2D 527 541 528 542 ! 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 543 DO_2D_00_00 544 zbt = zbet(ji,jj-1) 545 zbt1 = 1.0 - zbet(ji,jj-1) 546 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 547 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 548 zalf1 = 1.0 - zalf 549 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 550 ! 551 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 552 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 553 & + zbt1 * psy(ji,jj,jl) 554 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 555 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 556 & + zbt1 * psyy(ji,jj,jl) 557 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 558 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 559 & + zbt1 * psxy(ji,jj,jl) 560 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 561 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 562 END_2D 563 564 DO_2D_00_00 565 zbt = zbet(ji,jj) 566 zbt1 = 1.0 - zbet(ji,jj) 567 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 568 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 569 zalf1 = 1.0 - zalf 570 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 571 ! 572 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 573 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 574 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 575 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 576 & + ( zalf1 - zalf ) * ztemp ) ) 577 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 578 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 579 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 580 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 581 END_2D 572 582 573 583 END DO … … 579 589 ! 580 590 END SUBROUTINE adv_y 591 592 593 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 594 !!------------------------------------------------------------------- 595 !! *** ROUTINE Hbig *** 596 !! 597 !! ** Purpose : Thickness correction in case advection scheme creates 598 !! abnormally tick ice or snow 599 !! 600 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points 601 !! (before advection) and reduce it by adapting ice concentration 602 !! 2- check whether snow thickness is larger than the surrounding 9-points 603 !! (before advection) and reduce it by sending the excess in the ocean 604 !! 605 !! ** input : Max thickness of the surrounding 9-points 606 !!------------------------------------------------------------------- 607 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 608 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 609 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip 610 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 611 ! 612 INTEGER :: ji, jj, jl ! dummy loop indices 613 REAL(wp) :: z1_dt, zhip, zhi, zhs, zfra 614 !!------------------------------------------------------------------- 615 ! 616 z1_dt = 1._wp / pdt 617 ! 618 DO jl = 1, jpl 619 620 DO_2D_11_11 621 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 622 ! 623 ! ! -- check h_ip -- ! 624 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 625 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 626 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 627 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 628 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 629 ENDIF 630 ENDIF 631 ! 632 ! ! -- check h_i -- ! 633 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 634 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 635 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 636 pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 637 ENDIF 638 ! 639 ! ! -- check h_s -- ! 640 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 641 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 642 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 643 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 644 ! 645 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 646 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 647 ! 648 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 649 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 650 ENDIF 651 ! 652 ENDIF 653 END_2D 654 END DO 655 ! 656 END SUBROUTINE Hbig 581 657 582 658 … … 608 684 ! -- check snow load -- ! 609 685 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 ! 686 DO_2D_11_11 687 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 688 ! 689 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos ) 690 ! 691 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 692 ! put snow excess in the ocean 693 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 694 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 695 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 696 ! correct snow volume and heat content 697 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 698 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 626 699 ENDIF 627 END DO 628 END DO 700 ! 701 ENDIF 702 END_2D 629 703 END DO 630 704 ! … … 645 719 ! 646 720 ! !* 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) , & 721 ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 649 722 & sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , & 650 723 & sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , & … … 692 765 ! !==========================! 693 766 ! 694 IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sx opw' , ldstop = .FALSE. ) ! file exist: id1>0767 IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. ) ! file exist: id1>0 695 768 ELSE ; id1 = 0 ! no restart: id1=0 696 769 ENDIF … … 728 801 CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 729 802 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 803 ! ! snow layers heat content 737 804 DO jk = 1, nlay_s … … 776 843 sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity 777 844 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 845 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 780 846 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content … … 825 891 CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 826 892 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 893 ! ! snow layers heat content 834 894 DO jk = 1, nlay_s
Note: See TracChangeset
for help on using the changeset viewer.
