- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/icedyn_adv_pra.F90
r10425 r13463 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 20 21 USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call 22 USE icevar ! sea-ice: operations 21 23 ! 22 24 USE in_out_manager ! I/O manager … … 25 27 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 26 28 USE lbclnk ! lateral boundary conditions (or mpp links) 27 USE prtctl ! Print control28 29 29 30 IMPLICIT NONE … … 36 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxice, syice, sxxice, syyice, sxyice ! ice thickness 37 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsn , sysn , sxxsn , syysn , sxysn ! snow thickness 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! lead fraction39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! ice concentration 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 :: 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 86 INTEGER :: ji,jj, jk, jl, jt ! dummy loop indices 87 INTEGER :: icycle ! number of sub-timestep for the advection 88 REAL(wp) :: zdt ! - - 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 93 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 94 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 95 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp 96 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 97 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei 92 98 !!---------------------------------------------------------------------- 93 99 ! 94 100 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 95 101 ! 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 ) 102 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 103 DO jl = 1, jpl 104 DO_2D( 0, 0, 0, 0 ) 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.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 120 ! 121 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! 122 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 123 ! this should not affect too much the stability 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(:,:) ) ) 104 126 105 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 106 ELSE ; initad = 1 ; zusnit = 1.0_wp 127 ! non-blocking global communication send zcflnow and receive zcflprv 128 CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 129 130 IF( zcflprv(1) > .5 ) THEN ; icycle = 2 131 ELSE ; icycle = 1 107 132 ENDIF 133 zdt = rDt_ice / REAL(icycle) 108 134 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 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 154 END DO 155 DO jk = 1, nlay_i 156 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 157 END DO 158 IF ( ln_pnd_H12 ) THEN 159 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 160 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 161 ENDIF 122 162 END DO 123 DO jk = 1, nlay_i 124 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 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,:) ) 183 END DO 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,:) ) 221 END DO 222 IF ( ln_pnd_H12 ) THEN 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) 247 ENDIF 125 248 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 129 ENDIF 249 ! 250 ! derive open water from ice concentration 251 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 252 DO_2D( 0, 0, 0, 0 ) 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.0_wp ) 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 ! 130 270 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 271 ! 264 272 IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file … … 267 275 268 276 269 SUBROUTINE adv_x( pd f, put , pcrh, psm , ps0 , &277 SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , & 270 278 & psx, psxx, psy , psyy, psxy ) 271 279 !!---------------------------------------------------------------------- … … 275 283 !! variable on x axis 276 284 !!---------------------------------------------------------------------- 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 moments285 REAL(wp) , INTENT(in ) :: pdt ! the time step 286 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 287 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 288 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 289 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 290 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 291 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 284 292 !! 285 INTEGER :: ji, jj 286 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! local scalars293 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 294 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 287 295 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 288 296 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 291 299 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 292 300 !----------------------------------------------------------------------- 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 ! 302 jcat = SIZE( ps0 , 3 ) ! size of input arrays 303 ! 304 DO jl = 1, jcat ! loop on categories 305 ! 306 ! Limitation of moments. 307 DO_2D( 0, 0, 1, 1 ) 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) ) 301 312 zs1max = 1.5 * zslpmax 302 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj ) ) )313 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 303 314 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 304 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj ) ) )315 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 305 316 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 306 317 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 314 END DO 315 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+1 320 DO jj = 1, jpj ! Flux from i to i+1 WHEN u GT 0 321 DO ji = 1, jpi 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 325 326 ! Calculate fluxes and moments between boxes i<-->i+1 327 DO_2D( 0, 0, 1, 1 ) 322 328 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)329 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 324 330 zalfq = zalf * zalf 325 331 zalf1 = 1.0 - zalf 326 332 zalf1q = zalf1 * zalf1 327 333 ! 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 )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) 335 341 336 342 ! 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 DO 345 END DO 346 347 DO jj = 1, jpjm1 ! Flux from i+1 to i when u LT 0. 348 DO ji = 1, fs_jpim1 349 zalf = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 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( 0, 0, 1, 0 ) 353 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 350 354 zalg (ji,jj) = zalf 351 355 zalfq = zalf * zalf … … 355 359 zalg1q(ji,jj) = zalf1q 356 360 ! 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) * zalfq 361 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 DO 365 END DO 366 367 DO jj = 2, jpjm1 ! Readjust moments remaining in the box. 368 DO ji = fs_2, fs_jpim1 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( 0, 0, 0, 0 ) 369 372 zbt = zbet(ji-1,jj) 370 373 zbt1 = 1.0 - zbet(ji-1,jj) 371 374 ! 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 DO 380 END DO 381 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_jpim1 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 383 384 ! Put the temporary moments into appropriate neighboring boxes. 385 DO_2D( 0, 0, 0, 0 ) 385 386 zbt = zbet(ji-1,jj) 386 387 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 - zalf 390 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 DO 403 END DO 404 405 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0. 406 DO ji = fs_2, fs_jpim1 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( 0, 0, 0, 0 ) 407 406 zbt = zbet(ji,jj) 408 407 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 - zalf 412 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 DO 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 423 424 424 END DO 425 425 426 426 !-- 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 427 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 428 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 429 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 438 430 ! 439 431 END SUBROUTINE adv_x 440 432 441 433 442 SUBROUTINE adv_y( pd f, pvt , pcrh, psm , ps0 , &434 SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , & 443 435 & psx, psxx, psy , psyy, psxy ) 444 436 !!--------------------------------------------------------------------- … … 448 440 !! variable on y axis 449 441 !!--------------------------------------------------------------------- 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 moments457 !! 458 INTEGER :: ji, jj 459 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! temporary scalars442 REAL(wp) , INTENT(in ) :: pdt ! time step 443 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 444 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 445 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 446 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 447 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 448 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 449 !! 450 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 451 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 460 452 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 461 453 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 464 456 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 465 457 !--------------------------------------------------------------------- 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) ) 458 ! 459 jcat = SIZE( ps0 , 3 ) ! size of input arrays 460 ! 461 DO jl = 1, jcat ! loop on categories 462 ! 463 ! Limitation of moments. 464 DO_2D( 1, 1, 0, 0 ) 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) ) 474 469 zs1max = 1.5 * zslpmax 475 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj ) ) )470 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 476 471 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 477 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj ) ) )472 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 478 473 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 479 474 ! 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 487 END DO 488 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 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 482 483 ! Calculate fluxes and moments between boxes j<-->j+1 484 DO_2D( 1, 1, 0, 0 ) 495 485 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)486 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 497 487 zalfq = zalf * zalf 498 488 zalf1 = 1.0 - zalf 499 489 zalf1q = zalf1 * zalf1 500 490 ! 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 )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) 508 498 ! 509 499 ! Readjust moments remaining in the box. 510 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj) 511 ps0 (ji,jj) = ps0 (ji,jj) - zf0(ji,jj) 512 psy (ji,jj) = zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 513 psyy(ji,jj) = zalf1 * zalf1q * psyy(ji,jj) 514 psx (ji,jj) = psx (ji,jj) - zfx(ji,jj) 515 psxx(ji,jj) = psxx(ji,jj) - zfxx(ji,jj) 516 psxy(ji,jj) = zalf1q * psxy(ji,jj) 517 END DO 518 END DO 519 ! 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) 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( 1, 0, 0, 0 ) 510 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 523 511 zalg (ji,jj) = zalf 524 512 zalfq = zalf * zalf … … 528 516 zalg1q(ji,jj) = zalf1q 529 517 ! 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 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 527 528 ! Readjust moments remaining in the box. 529 DO_2D( 0, 0, 0, 0 ) 543 530 zbt = zbet(ji,jj-1) 544 531 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 545 532 ! 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 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 541 542 ! Put the temporary moments into appropriate neighboring boxes. 543 DO_2D( 0, 0, 0, 0 ) 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( 0, 0, 0, 0 ) 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 582 554 583 END DO 555 584 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 585 !-- Lateral boundary conditions 586 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 587 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 588 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 589 ! 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( 1, 1, 1, 1 ) 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 657 658 659 SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 660 !!------------------------------------------------------------------- 661 !! *** ROUTINE Hsnow *** 662 !! 663 !! ** Purpose : 1- Check snow load after advection 664 !! 2- Correct pond concentration to avoid a_ip > a_i 665 !! 666 !! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface 667 !! then put the snow excess in the ocean 668 !! 669 !! ** Notes : This correction is crucial because of the call to routine icecor afterwards 670 !! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 671 !! make the snow very thick (if concentration decreases drastically) 672 !! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 673 !!------------------------------------------------------------------- 674 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 675 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip 676 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 677 ! 678 INTEGER :: ji, jj, jl ! dummy loop indices 679 REAL(wp) :: z1_dt, zvs_excess, zfra 680 !!------------------------------------------------------------------- 681 ! 682 z1_dt = 1._wp / pdt 683 ! 684 ! -- check snow load -- ! 685 DO jl = 1, jpl 686 DO_2D( 1, 1, 1, 1 ) 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 699 ENDIF 700 ! 701 ENDIF 702 END_2D 578 703 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 704 ! 705 !-- correct pond concentration to avoid a_ip > a_i -- ! 706 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 707 ! 708 END SUBROUTINE Hsnow 614 709 615 710 … … 624 719 ! 625 720 ! !* 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) , & 721 ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 628 722 & sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , & 629 723 & sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , & … … 652 746 !! *** ROUTINE adv_pra_rst *** 653 747 !! 654 !! ** Purpose : Read or write RHGfile in restart file748 !! ** Purpose : Read or write file in restart file 655 749 !! 656 750 !! ** Method : use of IOM library … … 671 765 ! !==========================! 672 766 ! 673 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 674 768 ELSE ; id1 = 0 ! no restart: id1=0 675 769 ENDIF … … 678 772 ! 679 773 ! ! ice thickness 680 CALL iom_get( numrir, jpdom_auto glo, 'sxice' , sxice )681 CALL iom_get( numrir, jpdom_auto glo, 'syice' , syice )682 CALL iom_get( numrir, jpdom_auto glo, 'sxxice', sxxice )683 CALL iom_get( numrir, jpdom_auto glo, 'syyice', syyice )684 CALL iom_get( numrir, jpdom_auto glo, 'sxyice', sxyice )774 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice ) 775 CALL iom_get( numrir, jpdom_auto, 'syice' , syice ) 776 CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 777 CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 778 CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 685 779 ! ! snow thickness 686 CALL iom_get( numrir, jpdom_auto glo, 'sxsn' , sxsn )687 CALL iom_get( numrir, jpdom_auto glo, 'sysn' , sysn )688 CALL iom_get( numrir, jpdom_auto glo, 'sxxsn' , sxxsn )689 CALL iom_get( numrir, jpdom_auto glo, 'syysn' , syysn )690 CALL iom_get( numrir, jpdom_auto glo, 'sxysn' , sxysn )691 ! ! lead fraction692 CALL iom_get( numrir, jpdom_auto glo, 'sxa' , sxa )693 CALL iom_get( numrir, jpdom_auto glo, 'sya' , sya )694 CALL iom_get( numrir, jpdom_auto glo, 'sxxa' , sxxa )695 CALL iom_get( numrir, jpdom_auto glo, 'syya' , syya )696 CALL iom_get( numrir, jpdom_auto glo, 'sxya' , sxya )780 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn ) 781 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn ) 782 CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn ) 783 CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn ) 784 CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn ) 785 ! ! ice concentration 786 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa ) 787 CALL iom_get( numrir, jpdom_auto, 'sya' , sya ) 788 CALL iom_get( numrir, jpdom_auto, 'sxxa' , sxxa ) 789 CALL iom_get( numrir, jpdom_auto, 'syya' , syya ) 790 CALL iom_get( numrir, jpdom_auto, 'sxya' , sxya ) 697 791 ! ! ice salinity 698 CALL iom_get( numrir, jpdom_auto glo, 'sxsal' , sxsal )699 CALL iom_get( numrir, jpdom_auto glo, 'sysal' , sysal )700 CALL iom_get( numrir, jpdom_auto glo, 'sxxsal', sxxsal )701 CALL iom_get( numrir, jpdom_auto glo, 'syysal', syysal )702 CALL iom_get( numrir, jpdom_auto glo, 'sxysal', sxysal )792 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal ) 793 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal ) 794 CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 795 CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 796 CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 703 797 ! ! ice age 704 CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage ) 705 CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage ) 706 CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage ) 707 CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 708 CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage ) 709 ! ! open water in sea ice 710 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 ) 798 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage ) 799 CALL iom_get( numrir, jpdom_auto, 'syage' , syage ) 800 CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 801 CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) 802 CALL iom_get( numrir, jpdom_auto, 'sxyage', sxyage ) 715 803 ! ! snow layers heat content 716 804 DO jk = 1, nlay_s 717 805 WRITE(zchar1,'(I2.2)') jk 718 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)719 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:)720 znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:)721 znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:)722 znam = 'sxyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:)806 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:) 807 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:) 808 znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:) 809 znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:) 810 znam = 'sxyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:) 723 811 END DO 724 812 ! ! ice layers heat content 725 813 DO jk = 1, nlay_i 726 814 WRITE(zchar1,'(I2.2)') jk 727 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:)728 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:)729 znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:)730 znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:)731 znam = 'sxye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxye(:,:,jk,:) = z3d(:,:,:)815 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:) 816 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:) 817 znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:) 818 znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:) 819 znam = 'sxye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxye(:,:,jk,:) = z3d(:,:,:) 732 820 END DO 733 821 ! 734 822 IF( ln_pnd_H12 ) THEN ! melt pond fraction 735 CALL iom_get( numrir, jpdom_auto glo, 'sxap' , sxap )736 CALL iom_get( numrir, jpdom_auto glo, 'syap' , syap )737 CALL iom_get( numrir, jpdom_auto glo, 'sxxap', sxxap )738 CALL iom_get( numrir, jpdom_auto glo, 'syyap', syyap )739 CALL iom_get( numrir, jpdom_auto glo, 'sxyap', sxyap )823 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap ) 824 CALL iom_get( numrir, jpdom_auto, 'syap' , syap ) 825 CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 826 CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 827 CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 740 828 ! ! melt pond volume 741 CALL iom_get( numrir, jpdom_auto glo, 'sxvp' , sxvp )742 CALL iom_get( numrir, jpdom_auto glo, 'syvp' , syvp )743 CALL iom_get( numrir, jpdom_auto glo, 'sxxvp', sxxvp )744 CALL iom_get( numrir, jpdom_auto glo, 'syyvp', syyvp )745 CALL iom_get( numrir, jpdom_auto glo, 'sxyvp', sxyvp )829 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp ) 830 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp ) 831 CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 832 CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 833 CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 746 834 ENDIF 747 835 ! … … 752 840 sxice = 0._wp ; syice = 0._wp ; sxxice = 0._wp ; syyice = 0._wp ; sxyice = 0._wp ! ice thickness 753 841 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 fraction842 sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! ice concentration 755 843 sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity 756 844 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 845 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 759 846 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content … … 786 873 CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn ) 787 874 CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn ) 788 ! ! lead fraction875 ! ! ice concentration 789 876 CALL iom_rstput( iter, nitrst, numriw, 'sxa' , sxa ) 790 877 CALL iom_rstput( iter, nitrst, numriw, 'sya' , sya ) … … 804 891 CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 805 892 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 893 ! ! snow layers heat content 813 894 DO jk = 1, nlay_s
Note: See TracChangeset
for help on using the changeset viewer.