Changeset 5123 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
- Timestamp:
- 2015-03-04T17:06:03+01:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4990 r5123 17 17 USE dom_oce ! ocean domain 18 18 USE sbc_oce ! ocean surface boundary condition 19 USE par_ice ! ice parameter20 19 USE dom_ice ! ice domain 21 20 USE ice ! ice variables 22 21 USE limadv ! ice advection 23 22 USE limhdf ! ice horizontal diffusion 23 USE limvar ! 24 ! 24 25 USE in_out_manager ! I/O manager 25 26 USE lbclnk ! lateral boundary conditions -- MPP exchanges … … 28 29 USE prtctl ! Print control 29 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 USE limvar ! clem for ice thickness correction 31 USE timing ! Timing 31 USE timing ! Timing 32 32 USE limcons ! conservation tests 33 USE limctl ! control prints 33 34 34 35 IMPLICIT NONE 35 36 PRIVATE 36 37 37 PUBLIC lim_trp ! called by ice_step 38 PUBLIC lim_trp ! called by sbcice_lim 39 40 INTEGER :: ncfl ! number of ice time step with CFL>1/2 38 41 39 42 !! * Substitution … … 58 61 !! ** action : 59 62 !!--------------------------------------------------------------------- 60 INTEGER, INTENT(in) :: kt ! number of iteration63 INTEGER, INTENT(in) :: kt ! number of iteration 61 64 ! 62 INTEGER :: ji, jj, jk, jl, j n! dummy loop indices65 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 63 66 INTEGER :: initad ! number of sub-timestep for the advection 64 67 REAL(wp) :: zcfl , zusnit ! - - 68 CHARACTER(len=80) :: cltmp 65 69 ! 66 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 70 REAL(wp), POINTER, DIMENSION(:,:) :: zsm, zs0at 71 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ow 68 73 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 69 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 70 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 71 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 72 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 73 ! 74 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 75 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhimax ! old ice thickness 76 REAL(wp), POINTER, DIMENSION(:,:) :: zatold, zeiold, zesold ! old concentration, enthalpies 77 REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei 78 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 75 79 !!--------------------------------------------------------------------- 76 80 IF( nn_timing == 1 ) CALL timing_start('limtrp') 77 81 78 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )79 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi)80 CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e)81 82 CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem82 CALL wrk_alloc( jpi,jpj, zsm, zs0at, zatold, zeiold, zesold ) 83 CALL wrk_alloc( jpi,jpj,jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e ) 84 CALL wrk_alloc( jpi,jpj,1, zs0ow ) 85 CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 86 CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold, zsmvold ) 83 87 84 88 IF( numit == nstart .AND. lwp ) THEN … … 88 92 ENDIF 89 93 WRITE(numout,*) '~~~~~~~~~~~~' 94 ncfl = 0 ! nb of time step with CFL > 1/2 90 95 ENDIF 96 97 zsm(:,:) = e12t(:,:) 91 98 92 zsm(:,:) = area(:,:)93 94 99 ! !-------------------------------------! 95 100 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! … … 97 102 98 103 ! conservation test 99 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)100 101 ! mass and salt flux init (clem)104 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 105 106 ! mass and salt flux init 102 107 zviold(:,:,:) = v_i(:,:,:) 103 zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 104 zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 105 106 !--- Thickness correction init. (clem) ------------------------------- 108 zvsold(:,:,:) = v_s(:,:,:) 109 zsmvold(:,:,:) = smv_i(:,:,:) 110 zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 111 zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 112 113 !--- Thickness correction init. ------------------------------- 107 114 CALL lim_var_glo2eqv 108 za iold(:,:,:) = a_i(:,:,:)115 zatold(:,:) = SUM( a_i(:,:,:), dim=3 ) 109 116 !--------------------------------------------------------------------- 110 117 ! Record max of the surrounding ice thicknesses for correction in limupdate … … 116 123 DO ji = 2, jpim1 117 124 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) ) 118 !zhimax(ji,jj,jl) = ( ht_i(ji ,jj ,jl) * tmask(ji, jj ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) &119 ! & + ht_i(ji-1,jj ,jl) * tmask(ji-1,jj ,1) + ht_i(ji ,jj-1,jl) * tmask(ji ,jj-1,1) &120 ! & + ht_i(ji+1,jj ,jl) * tmask(ji+1,jj ,1) + ht_i(ji ,jj+1,jl) * tmask(ji ,jj+1,1) &121 ! & + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) )122 125 END DO 123 126 END DO … … 125 128 END DO 126 129 130 !=============================! 131 !== Prather scheme ==! 132 !=============================! 133 134 ! If ice drift field is too fast, use an appropriate time step for advection. 135 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability 136 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 137 IF(lk_mpp ) CALL mpp_max( zcfl ) 138 139 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 140 ELSE ; initad = 1 ; zusnit = 1.0_wp 141 ENDIF 142 143 IF( zcfl > 0.5_wp .AND. lwp ) ncfl = ncfl + 1 144 IF( numit == nlast .AND. lwp ) THEN 145 IF( ncfl > 0 ) THEN 146 WRITE(cltmp,'(i6.1)') ncfl 147 CALL ctl_stop('STOP',TRIM(cltmp) ) 148 CALL ctl_warn( 'lim_trp: ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 149 ELSE 150 WRITE(numout,*) 'lim_trp : CFL criteria for ice advection is always smaller than 1/2 ' 151 ENDIF 152 ENDIF 153 127 154 !------------------------- 128 155 ! transported fields 129 156 !------------------------- 130 ! Snow vol, ice vol, salt and age contents, area 131 zs0ow(:,:) = ato_i(:,:) * area(:,:) ! Open water area 132 DO jl = 1, jpl 133 zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume 134 zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume 135 zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area 136 zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content 137 zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content 138 zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content 139 zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content 140 END DO 141 142 !-------------------------- 143 ! Advection of Ice fields (Prather scheme) 144 !-------------------------- 145 ! If ice drift field is too fast, use an appropriate time step for advection. 146 ! CFL test for stability 147 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) ) 148 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 149 IF(lk_mpp ) CALL mpp_max( zcfl ) 150 !!gm more readability: 151 ! IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 152 ! ELSE ; initad = 1 ; zusnit = 1.0_wp 153 ! ENDIF 154 !!gm end 155 initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) ) 156 zusnit = 1.0 / REAL( initad ) 157 IF( zcfl > 0.5 .AND. lwp ) & 158 WRITE(numout,*) 'lim_trp : CFL violation at day ', nday, ', cfl = ', zcfl, & 159 & ': the ice time stepping is split in two' 157 zs0ow(:,:,1) = ato_i(:,:) * e12t(:,:) ! Open water area 158 DO jl = 1, jpl 159 zs0sn (:,:,jl) = v_s (:,:,jl) * e12t(:,:) ! Snow volume 160 zs0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume 161 zs0a (:,:,jl) = a_i (:,:,jl) * e12t(:,:) ! Ice area 162 zs0sm (:,:,jl) = smv_i(:,:,jl) * e12t(:,:) ! Salt content 163 zs0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content 164 zs0c0 (:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content 165 DO jk = 1, nlay_i 166 zs0e (:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content 167 END DO 168 END DO 169 160 170 161 171 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 162 DO j n = 1,initad163 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area164 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )165 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,: ), sxopw(:,:), &166 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )172 DO jt = 1, initad 173 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:), & !--- ice open water area 174 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 175 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:), & 176 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 167 177 DO jl = 1, jpl 168 CALL lim_adv_x( zusnit, u_ice, 1._wp 178 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 169 179 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 170 180 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 171 181 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 172 CALL lim_adv_x( zusnit, u_ice, 1._wp 182 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 173 183 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 174 184 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 175 185 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 176 CALL lim_adv_x( zusnit, u_ice, 1._wp 186 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 177 187 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 178 188 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 179 189 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 180 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), &!--- ice age ---190 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 181 191 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 182 192 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 183 193 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 184 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), &!--- ice concentrations ---194 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 185 195 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 186 196 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 187 197 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 188 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &!--- snow heat contents ---198 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 189 199 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 190 200 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 191 201 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 192 DO jk = 1, nlay_i !--- ice heat contents ---193 CALL lim_adv_x( zusnit, u_ice, 1._wp 202 DO jk = 1, nlay_i !--- ice heat contents --- 203 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 194 204 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 195 205 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) … … 201 211 END DO 202 212 ELSE 203 DO j n= 1, initad204 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area205 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )206 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,: ), sxopw(:,:), &207 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )213 DO jt = 1, initad 214 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:), & !--- ice open water area 215 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 216 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:), & 217 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 208 218 DO jl = 1, jpl 209 CALL lim_adv_y( zusnit, v_ice, 1._wp 219 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 210 220 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 211 221 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 212 222 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 213 CALL lim_adv_y( zusnit, v_ice, 1._wp 223 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 214 224 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 215 225 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 216 226 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 217 CALL lim_adv_y( zusnit, v_ice, 1._wp 227 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 218 228 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 219 229 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 220 230 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 221 231 222 CALL lim_adv_y( zusnit, v_ice, 1._wp 232 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 223 233 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 224 234 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 225 235 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 226 CALL lim_adv_y( zusnit, v_ice, 1._wp 236 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 227 237 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 228 238 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 229 239 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 230 CALL lim_adv_y( zusnit, v_ice, 1._wp 240 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 231 241 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 232 242 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 233 243 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 234 244 DO jk = 1, nlay_i !--- ice heat contents --- 235 CALL lim_adv_y( zusnit, v_ice, 1._wp 245 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 236 246 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 237 247 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) … … 247 257 ! Recover the properties from their contents 248 258 !------------------------------------------- 249 zs0ow(:,:) = zs0ow(:,:) / area(:,:) 250 DO jl = 1, jpl 251 zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:) 252 zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:) 253 zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:) 254 zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 255 zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:) 256 ! 259 ato_i(:,:) = zs0ow(:,:,1) * r1_e12t(:,:) 260 DO jl = 1, jpl 261 v_i (:,:,jl) = zs0ice(:,:,jl) * r1_e12t(:,:) 262 v_s (:,:,jl) = zs0sn (:,:,jl) * r1_e12t(:,:) 263 smv_i(:,:,jl) = zs0sm (:,:,jl) * r1_e12t(:,:) 264 oa_i (:,:,jl) = zs0oi (:,:,jl) * r1_e12t(:,:) 265 a_i (:,:,jl) = zs0a (:,:,jl) * r1_e12t(:,:) 266 e_s (:,:,1,jl) = zs0c0 (:,:,jl) * r1_e12t(:,:) 267 DO jk = 1, nlay_i 268 e_i (:,:,jk,jl) = zs0e (:,:,jk,jl) * r1_e12t(:,:) 269 END DO 270 END DO 271 272 at_i(:,:) = a_i(:,:,1) ! total ice fraction 273 DO jl = 2, jpl 274 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 257 275 END DO 258 276 259 277 !------------------------------------------------------------------------------! 260 ! 4)Diffusion of Ice fields278 ! Diffusion of Ice fields 261 279 !------------------------------------------------------------------------------! 262 280 281 ! 263 282 !-------------------------------- 264 283 ! diffusion of open water area 265 284 !-------------------------------- 266 zs0at(:,:) = zs0a(:,:,1) ! total ice fraction267 DO jl = 2, jpl268 zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)269 END DO270 !271 285 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points 272 286 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 273 287 DO ji = 1 , fs_jpim1 ! vector opt. 274 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0at(ji ,jj) ) ) ) &275 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)276 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0at(ji,jj ) ) ) ) &277 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)288 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji ,jj) ) ) ) & 289 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 290 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj ) ) ) ) & 291 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 278 292 END DO 279 293 END DO 280 294 ! 281 CALL lim_hdf( zs0ow(:,:) ) ! Diffusion295 CALL lim_hdf( ato_i (:,:) ) ! Diffusion 282 296 283 297 !------------------------------------ … … 288 302 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 289 303 DO ji = 1 , fs_jpim1 ! vector opt. 290 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0a(ji ,jj,jl) ) ) ) &291 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)292 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0a(ji,jj ,jl) ) ) ) &293 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)294 END DO 295 END DO 296 297 CALL lim_hdf( zs0ice (:,:,jl) )298 CALL lim_hdf( zs0sn (:,:,jl) )299 CALL lim_hdf( zs0sm (:,:,jl) )300 CALL lim_hdf( zs0oi (:,:,jl) )301 CALL lim_hdf( zs0a (:,:,jl) )302 CALL lim_hdf( zs0c0 (:,:,jl) )304 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji ,jj,jl) ) ) ) & 305 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 306 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj ,jl) ) ) ) & 307 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 308 END DO 309 END DO 310 311 CALL lim_hdf( v_i (:,:, jl) ) 312 CALL lim_hdf( v_s (:,:, jl) ) 313 CALL lim_hdf( smv_i(:,:, jl) ) 314 CALL lim_hdf( oa_i (:,:, jl) ) 315 CALL lim_hdf( a_i (:,:, jl) ) 316 CALL lim_hdf( e_s (:,:,1,jl) ) 303 317 DO jk = 1, nlay_i 304 CALL lim_hdf( zs0e(:,:,jk,jl) )318 CALL lim_hdf( e_i(:,:,jk,jl) ) 305 319 END DO 306 320 END DO 307 321 308 322 !------------------------------------------------------------------------------! 309 ! 5) Update andlimit ice properties after transport323 ! limit ice properties after transport 310 324 !------------------------------------------------------------------------------! 311 312 !-------------------------------------------------- 313 ! 5.1) Recover mean values over the grid squares. 314 !-------------------------------------------------- 315 zs0at(:,:) = 0._wp 325 !!gm & cr : MAX should not be active if adv scheme is positive ! 316 326 DO jl = 1, jpl 317 327 DO jj = 1, jpj 318 328 DO ji = 1, jpi 319 zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) ) 320 zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) ) 321 zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) ) 322 zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) ) 323 zs0a (ji,jj,jl) = MAX( 0._wp, zs0a (ji,jj,jl) ) 324 zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) ) 325 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 326 END DO 327 END DO 328 END DO 329 330 !--------------------------------------------------------- 331 ! 5.2) Update and mask variables 332 !--------------------------------------------------------- 333 DO jl = 1, jpl 334 DO jj = 1, jpj 335 DO ji = 1, jpi 336 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 337 338 zvi = zs0ice(ji,jj,jl) 339 zvs = zs0sn (ji,jj,jl) 340 zes = zs0c0 (ji,jj,jl) 341 zsmv = zs0sm (ji,jj,jl) 342 ! 343 ! Remove very small areas 344 v_s(ji,jj,jl) = rswitch * zs0sn (ji,jj,jl) 345 v_i(ji,jj,jl) = rswitch * zs0ice(ji,jj,jl) 346 a_i(ji,jj,jl) = rswitch * zs0a (ji,jj,jl) 347 e_s(ji,jj,1,jl) = rswitch * zs0c0 (ji,jj,jl) 348 ! Ice salinity and age 349 IF( num_sal == 2 ) THEN 350 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 351 ENDIF 352 oa_i(ji,jj,jl) = MAX( rswitch * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 353 354 ! Update fluxes 355 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 356 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 357 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 358 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 359 END DO 360 END DO 361 END DO 362 363 DO jl = 1, jpl 329 v_s (ji,jj,jl) = MAX( 0._wp, v_s (ji,jj,jl) ) 330 v_i (ji,jj,jl) = MAX( 0._wp, v_i (ji,jj,jl) ) 331 smv_i(ji,jj,jl) = MAX( 0._wp, smv_i(ji,jj,jl) ) 332 oa_i (ji,jj,jl) = MAX( 0._wp, oa_i (ji,jj,jl) ) 333 a_i (ji,jj,jl) = MAX( 0._wp, a_i (ji,jj,jl) ) 334 e_s (ji,jj,1,jl) = MAX( 0._wp, e_s (ji,jj,1,jl) ) 335 END DO 336 END DO 337 364 338 DO jk = 1, nlay_i 365 339 DO jj = 1, jpj 366 340 DO ji = 1, jpi 367 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10) )368 zei = zs0e(ji,jj,jk,jl)369 e_i(ji,jj,jk,jl) = rswitch * MAX( 0._wp, zs0e(ji,jj,jk,jl) )370 ! Update fluxes371 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0372 END DO !ji 373 END DO ! jj 374 END DO ! jk375 END DO ! jl376 377 !--- Thickness correction in case too high (clem)--------------------------------------------------------341 e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) ) 342 END DO 343 END DO 344 END DO 345 END DO 346 !!gm & cr 347 348 ! zap small areas 349 CALL lim_var_zapsmall 350 351 !--- Thickness correction in case too high -------------------------------------------------------- 378 352 CALL lim_var_glo2eqv 379 353 DO jl = 1, jpl … … 388 362 zei = SUM( e_i(ji,jj,1:nlay_i,jl) ) 389 363 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 390 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)391 364 392 365 rswitch = 1._wp 393 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl)) < 0.80 ) .OR. &366 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 394 367 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 395 368 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) … … 413 386 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 414 387 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 415 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) *r1_rdtice ! W.m-2 <0416 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) *r1_rdtice ! W.m-2 <0388 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 389 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 417 390 ENDIF 391 418 392 END DO 419 393 END DO 420 394 END DO 421 395 ! ------------------------------------------------- 396 397 !-------------------------------------- 398 ! Impose a_i < amax in mono-category 399 !-------------------------------------- 400 ! 401 IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2 402 DO jj = 1, jpj 403 DO ji = 1, jpi 404 a_i(ji,jj,1) = MIN( a_i(ji,jj,1), rn_amax ) 405 END DO 406 END DO 407 ENDIF 422 408 423 409 ! --- diags --- 424 410 DO jj = 1, jpj 425 411 DO ji = 1, jpi 426 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 427 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 428 429 diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 430 diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 412 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 413 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 414 415 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 416 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 417 diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 431 418 END DO 432 419 END DO … … 454 441 ! open water = 1 if at_i=0 455 442 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 456 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj)443 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 457 444 END DO 458 445 END DO … … 463 450 ENDIF 464 451 465 IF(ln_ctl) THEN ! Control print 466 CALL prt_ctl_info(' ') 467 CALL prt_ctl_info(' - Cell values : ') 468 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 469 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp : cell area :') 470 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp : at_i :') 471 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp : vt_i :') 472 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp : vt_s :') 473 DO jl = 1, jpl 474 CALL prt_ctl_info(' ') 475 CALL prt_ctl_info(' - Category : ', ivar1=jl) 476 CALL prt_ctl_info(' ~~~~~~~~~~') 477 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_trp : a_i : ') 478 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_trp : ht_i : ') 479 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_trp : ht_s : ') 480 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_trp : v_i : ') 481 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_trp : v_s : ') 482 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_trp : e_s : ') 483 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_trp : t_su : ') 484 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_trp : t_snow : ') 485 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_trp : sm_i : ') 486 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_trp : smv_i : ') 487 DO jk = 1, nlay_i 488 CALL prt_ctl_info(' ') 489 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 490 CALL prt_ctl_info(' ~~~~~~~') 491 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp : t_i : ') 492 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp : e_i : ') 493 END DO 494 END DO 495 ENDIF 452 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting 453 454 ! ------------------------------------------------- 455 ! control prints 456 ! ------------------------------------------------- 457 IF( ln_nicep ) CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' ) 496 458 ! 497 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )498 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi)499 CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e)500 501 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax ) ! clem459 CALL wrk_dealloc( jpi,jpj, zsm, zs0at, zatold, zeiold, zesold ) 460 CALL wrk_dealloc( jpi,jpj,jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e ) 461 CALL wrk_dealloc( jpi,jpj,1, zs0ow ) 462 CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 463 CALL wrk_dealloc( jpi,jpj,jpl, zviold, zvsold, zhimax, zsmvold ) 502 464 ! 503 465 IF( nn_timing == 1 ) CALL timing_stop('limtrp') 466 504 467 END SUBROUTINE lim_trp 505 468 … … 512 475 END SUBROUTINE lim_trp 513 476 #endif 514 515 477 !!====================================================================== 516 478 END MODULE limtrp
Note: See TracChangeset
for help on using the changeset viewer.
