Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
- Timestamp:
- 2015-05-12T12:37:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4333 r5260 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 USE limcons ! conservation tests 33 USE limctl ! control prints 32 34 33 35 IMPLICIT NONE 34 36 PRIVATE 35 37 36 PUBLIC lim_trp ! called by ice_step 37 38 REAL(wp) :: epsi10 = 1.e-10_wp 39 REAL(wp) :: rzero = 0._wp 40 REAL(wp) :: rone = 1._wp 38 PUBLIC lim_trp ! called by sbcice_lim 39 40 INTEGER :: ncfl ! number of ice time step with CFL>1/2 41 41 42 42 !! * Substitution … … 61 61 !! ** action : 62 62 !!--------------------------------------------------------------------- 63 INTEGER, INTENT(in) :: kt ! number of iteration63 INTEGER, INTENT(in) :: kt ! number of iteration 64 64 ! 65 INTEGER :: ji, jj, jk, jl, layer! dummy loop indices65 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 66 66 INTEGER :: initad ! number of sub-timestep for the advection 67 INTEGER :: ierr ! error status 68 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 69 REAL(wp) :: zusvosn, zusvoic, zbigval ! - - 70 REAL(wp) :: zcfl , zusnit ! - - 71 REAL(wp) :: ze , zsal , zage ! - - 67 REAL(wp) :: zcfl , zusnit ! - - 68 CHARACTER(len=80) :: cltmp 72 69 ! 73 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 75 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 76 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 77 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset) 78 ! mass and salt flux (clem) 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold ! old ice volume... 80 ! correct ice thickness (clem) 81 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 82 REAL(wp) :: zdv, zda, zvi, zvs, zsmv 70 REAL(wp), POINTER, DIMENSION(:,:) :: zsm 71 REAL(wp), POINTER, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0es , z0smi , z0oi 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: z0opw 73 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: z0ei 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 83 79 !!--------------------------------------------------------------------- 84 80 IF( nn_timing == 1 ) CALL timing_start('limtrp') 85 81 86 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow ) 87 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 88 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e ) 89 90 CALL wrk_alloc( jpi,jpj,jpl,zviold ) ! clem 91 CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 92 93 ! ------------------------------- 94 !- check conservation (C Rousset) 95 IF( ln_limdiahsb ) THEN 96 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 97 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 98 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 99 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 100 ENDIF 101 !- check conservation (C Rousset) 102 ! ------------------------------- 82 CALL wrk_alloc( jpi,jpj, zsm, zatold, zeiold, zesold ) 83 CALL wrk_alloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 84 CALL wrk_alloc( jpi,jpj,1, z0opw ) 85 CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 86 CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold, zsmvold ) 103 87 104 88 IF( numit == nstart .AND. lwp ) THEN … … 108 92 ENDIF 109 93 WRITE(numout,*) '~~~~~~~~~~~~' 94 ncfl = 0 ! nb of time step with CFL > 1/2 110 95 ENDIF 96 97 zsm(:,:) = e12t(:,:) 111 98 112 zsm(:,:) = area(:,:)113 114 99 ! !-------------------------------------! 115 100 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 116 101 ! !-------------------------------------! 117 ! mass and salt flux init (clem) 102 103 ! conservation test 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 118 107 zviold(:,:,:) = v_i(:,:,:) 119 120 !--- Thickness correction init. (clem) ------------------------------- 121 CALL lim_var_glo2eqv 122 zaiold(:,:,:) = a_i(:,:,:) 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. ------------------------------- 114 zatold(:,:) = SUM( a_i(:,:,:), dim=3 ) 115 DO jl = 1, jpl 116 DO jj = 1, jpj 117 DO ji = 1, jpi 118 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 119 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 120 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 121 END DO 122 END DO 123 END DO 123 124 !--------------------------------------------------------------------- 124 ! Record max of the surrounding ice thicknesses for correction in limupdate125 ! Record max of the surrounding ice thicknesses for correction 125 126 ! in case advection creates ice too thick. 126 127 !--------------------------------------------------------------------- 127 zhimax(:,:,:) = ht_i(:,:,:) 128 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 128 129 DO jl = 1, jpl 129 130 DO jj = 2, jpjm1 130 131 DO ji = 2, jpim1 131 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) ) 132 !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) & 133 ! & + ht_i(ji-1,jj ,jl) * tmask(ji-1,jj ,1) + ht_i(ji ,jj-1,jl) * tmask(ji ,jj-1,1) & 134 ! & + ht_i(ji+1,jj ,jl) * tmask(ji+1,jj ,1) + ht_i(ji ,jj+1,jl) * tmask(ji ,jj+1,1) & 135 ! & + 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) ) 132 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) 136 133 END DO 137 134 END DO … … 139 136 END DO 140 137 138 !=============================! 139 !== Prather scheme ==! 140 !=============================! 141 142 ! If ice drift field is too fast, use an appropriate time step for advection. 143 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability 144 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 145 IF(lk_mpp ) CALL mpp_max( zcfl ) 146 147 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 148 ELSE ; initad = 1 ; zusnit = 1.0_wp 149 ENDIF 150 151 IF( zcfl > 0.5_wp .AND. lwp ) ncfl = ncfl + 1 152 !! IF( lwp ) THEN 153 !! IF( ncfl > 0 ) THEN 154 !! WRITE(cltmp,'(i6.1)') ncfl 155 !! CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 156 !! ELSE 157 !! ! WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 158 !! ENDIF 159 !! ENDIF 160 141 161 !------------------------- 142 162 ! transported fields 143 163 !------------------------- 144 ! Snow vol, ice vol, salt and age contents, area 145 zs0ow(:,:) = ato_i(:,:) * area(:,:) ! Open water area 146 DO jl = 1, jpl 147 zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume 148 zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume 149 zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area 150 zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content 151 zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content 152 zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content 153 zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content 154 END DO 155 156 !-------------------------- 157 ! Advection of Ice fields (Prather scheme) 158 !-------------------------- 159 ! If ice drift field is too fast, use an appropriate time step for advection. 160 ! CFL test for stability 161 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) ) 162 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 163 IF(lk_mpp ) CALL mpp_max( zcfl ) 164 !!gm more readability: 165 ! IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 166 ! ELSE ; initad = 1 ; zusnit = 1.0_wp 167 ! ENDIF 168 !!gm end 169 initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 170 zusnit = 1.0 / REAL( initad ) 171 IF( zcfl > 0.5 .AND. lwp ) & 172 WRITE(numout,*) 'lim_trp : CFL violation at day ', nday, ', cfl = ', zcfl, & 173 & ': the ice time stepping is split in two' 164 z0opw(:,:,1) = ato_i(:,:) * e12t(:,:) ! Open water area 165 DO jl = 1, jpl 166 z0snw (:,:,jl) = v_s (:,:,jl) * e12t(:,:) ! Snow volume 167 z0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume 168 z0ai (:,:,jl) = a_i (:,:,jl) * e12t(:,:) ! Ice area 169 z0smi (:,:,jl) = smv_i(:,:,jl) * e12t(:,:) ! Salt content 170 z0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content 171 z0es (:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content 172 DO jk = 1, nlay_i 173 z0ei (:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content 174 END DO 175 END DO 176 174 177 175 178 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 176 DO j k = 1,initad177 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area178 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )179 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &180 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )179 DO jt = 1, initad 180 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 181 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 182 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), & 183 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 181 184 DO jl = 1, jpl 182 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---185 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 183 186 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 184 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &187 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & 185 188 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 186 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn(:,:,jl), sxsn (:,:,jl), & !--- snow volume ---189 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 187 190 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 188 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn(:,:,jl), sxsn (:,:,jl), &191 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & 189 192 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 190 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm(:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---193 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 191 194 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 192 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm(:,:,jl), sxsal(:,:,jl), &195 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & 193 196 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 194 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), &!--- ice age ---197 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 195 198 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi(:,:,jl), sxage(:,:,jl), &199 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & 197 200 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 198 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a (:,:,jl), sxa (:,:,jl), &!--- ice concentrations ---201 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 199 202 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 200 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a(:,:,jl), sxa (:,:,jl), &203 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & 201 204 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 202 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &!--- snow heat contents ---205 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 203 206 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 204 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0(:,:,jl), sxc0 (:,:,jl), &207 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & 205 208 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 206 DO layer = 1, nlay_i!--- ice heat contents ---207 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &208 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &209 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )210 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &211 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &212 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )209 DO jk = 1, nlay_i !--- ice heat contents --- 210 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 211 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 212 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 213 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 214 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 215 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 213 216 END DO 214 217 END DO 215 218 END DO 216 219 ELSE 217 DO j k= 1, initad218 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area219 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )220 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &221 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )220 DO jt = 1, initad 221 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 222 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 223 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), & 224 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 222 225 DO jl = 1, jpl 223 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---226 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 224 227 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 225 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &228 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & 226 229 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 227 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sn(:,:,jl), sxsn (:,:,jl), & !--- snow volume ---230 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 228 231 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 229 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn(:,:,jl), sxsn (:,:,jl), &232 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & 230 233 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 231 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sm(:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---234 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 232 235 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 233 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm(:,:,jl), sxsal(:,:,jl), &236 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & 234 237 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 235 236 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 238 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 237 239 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi(:,:,jl), sxage(:,:,jl), &240 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & 239 241 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 240 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0a(:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---242 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 241 243 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 242 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a(:,:,jl), sxa (:,:,jl), &244 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & 243 245 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 244 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0c0(:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---246 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 245 247 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 246 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0(:,:,jl), sxc0 (:,:,jl), &248 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & 247 249 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 248 DO layer= 1, nlay_i !--- ice heat contents ---249 CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &250 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &251 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )252 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &253 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &254 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )250 DO jk = 1, nlay_i !--- ice heat contents --- 251 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 252 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 253 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 254 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 255 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 256 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 255 257 END DO 256 258 END DO … … 261 263 ! Recover the properties from their contents 262 264 !------------------------------------------- 263 zs0ow(:,:) = zs0ow(:,:) / area(:,:)264 DO jl = 1, jpl 265 zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)266 zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)267 zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)268 zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)269 zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:)270 zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:)265 ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:) 266 DO jl = 1, jpl 267 v_i (:,:,jl) = z0ice(:,:,jl) * r1_e12t(:,:) 268 v_s (:,:,jl) = z0snw(:,:,jl) * r1_e12t(:,:) 269 smv_i(:,:,jl) = z0smi(:,:,jl) * r1_e12t(:,:) 270 oa_i (:,:,jl) = z0oi (:,:,jl) * r1_e12t(:,:) 271 a_i (:,:,jl) = z0ai (:,:,jl) * r1_e12t(:,:) 272 e_s (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:) 271 273 DO jk = 1, nlay_i 272 zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:) 273 END DO 274 e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:) 275 END DO 276 END DO 277 278 at_i(:,:) = a_i(:,:,1) ! total ice fraction 279 DO jl = 2, jpl 280 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 274 281 END DO 275 282 276 283 !------------------------------------------------------------------------------! 277 ! 4)Diffusion of Ice fields284 ! Diffusion of Ice fields 278 285 !------------------------------------------------------------------------------! 279 286 287 ! 280 288 !-------------------------------- 281 289 ! diffusion of open water area 282 290 !-------------------------------- 283 zs0at(:,:) = zs0a(:,:,1) ! total ice fraction284 DO jl = 2, jpl285 zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)286 END DO287 !288 291 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points 289 292 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 290 293 DO ji = 1 , fs_jpim1 ! vector opt. 291 pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji ,jj) ) ) ) &292 & * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)293 pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj ) ) ) ) &294 & * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)294 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji ,jj) ) ) ) & 295 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 296 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj ) ) ) ) & 297 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 295 298 END DO 296 299 END DO 297 300 ! 298 CALL lim_hdf( zs0ow (:,:) ) ! Diffusion301 CALL lim_hdf( ato_i (:,:) ) 299 302 300 303 !------------------------------------ … … 305 308 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 306 309 DO ji = 1 , fs_jpim1 ! vector opt. 307 pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji ,jj,jl) ) ) ) &308 & * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)309 pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj ,jl) ) ) ) &310 & * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)311 END DO 312 END DO 313 314 CALL lim_hdf( zs0ice (:,:,jl) )315 CALL lim_hdf( zs0sn (:,:,jl) )316 CALL lim_hdf( zs0sm (:,:,jl) )317 CALL lim_hdf( zs0oi (:,:,jl) )318 CALL lim_hdf( zs0a (:,:,jl) )319 CALL lim_hdf( zs0c0 (:,:,jl) )310 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji ,jj,jl) ) ) ) & 311 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 312 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj ,jl) ) ) ) & 313 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 314 END DO 315 END DO 316 317 CALL lim_hdf( v_i (:,:, jl) ) 318 CALL lim_hdf( v_s (:,:, jl) ) 319 CALL lim_hdf( smv_i(:,:, jl) ) 320 CALL lim_hdf( oa_i (:,:, jl) ) 321 CALL lim_hdf( a_i (:,:, jl) ) 322 CALL lim_hdf( e_s (:,:,1,jl) ) 320 323 DO jk = 1, nlay_i 321 CALL lim_hdf( zs0e(:,:,jk,jl) )324 CALL lim_hdf( e_i(:,:,jk,jl) ) 322 325 END DO 323 326 END DO 324 327 325 328 !------------------------------------------------------------------------------! 326 ! 5) Update andlimit ice properties after transport329 ! limit ice properties after transport 327 330 !------------------------------------------------------------------------------! 328 329 !-------------------------------------------------- 330 ! 5.1) Recover mean values over the grid squares. 331 !-------------------------------------------------- 332 zs0at(:,:) = 0._wp 331 !!gm & cr : MAX should not be active if adv scheme is positive ! 333 332 DO jl = 1, jpl 334 333 DO jj = 1, jpj 335 334 DO ji = 1, jpi 336 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) ) 337 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) ) 338 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) ) 339 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) ) 340 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) ) 341 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) ) 342 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 343 END DO 344 END DO 345 END DO 346 347 !--------------------------------------------------------- 348 ! 5.2) Snow thickness, Ice thickness, Ice concentrations 349 !--------------------------------------------------------- 350 DO jj = 1, jpj 351 DO ji = 1, jpi 352 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) ) 353 zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp ) 354 ato_i(ji,jj) = zs0ow(ji,jj) 355 END DO 356 END DO 357 358 DO jl = 1, jpl ! Remove very small areas 359 DO jj = 1, jpj 360 DO ji = 1, jpi 361 zvi = zs0ice(ji,jj,jl) 362 zvs = zs0sn(ji,jj,jl) 363 ! 364 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) ) 365 ! 366 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 367 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 368 ! 369 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 370 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 371 zindb = MAX( zindsn, zindic ) 372 ! 373 zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration 374 a_i (ji,jj,jl) = zs0a(ji,jj,jl) 375 v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 376 v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl) 377 ! 378 ! Update mass fluxes (clem) 379 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 380 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 381 END DO 382 END DO 383 END DO 384 385 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 386 CALL lim_var_glo2eqv 387 DO jl = 1, jpl 388 DO jj = 1, jpj 389 DO ji = 1, jpi 390 391 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 392 zvi = v_i(ji,jj,jl) 393 zvs = v_s(ji,jj,jl) 394 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 395 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 396 397 zindh = 1._wp 398 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 399 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 400 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 401 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 402 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 403 ELSE 404 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 405 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 406 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 407 ENDIF 408 409 ! small correction due to *zindh for a_i 410 v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) 411 v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl) 412 413 ! Update mass fluxes 414 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 415 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 416 417 ENDIF 418 419 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 420 421 END DO 422 END DO 423 END DO 424 425 ! --- 426 DO jj = 1, jpj 427 DO ji = 1, jpi 428 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless?? 429 END DO 430 END DO 431 432 !---------------------- 433 ! 5.3) Ice properties 434 !---------------------- 435 436 zbigval = 1.e+13 437 438 DO jl = 1, jpl 439 DO jj = 1, jpj 440 DO ji = 1, jpi 441 zsmv = zs0sm(ji,jj,jl) 442 443 ! Switches and dummy variables 444 zusvosn = 1.0/MAX( v_s(ji,jj,jl) , epsi10 ) 445 zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi10 ) 446 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 447 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 448 zindb = MAX( zindsn, zindic ) 449 450 ! Ice salinity and age 451 !clem zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj), zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl) 452 IF( num_sal == 2 ) THEN 453 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 454 ENDIF 455 456 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp ) * a_i(ji,jj,jl) 457 oa_i (ji,jj,jl) = zindic * zage 458 459 ! Snow heat content 460 ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval ) 461 e_s(ji,jj,1,jl) = zindsn * ze 462 463 ! Update salt fluxes (clem) 464 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 465 END DO !ji 466 END DO !jj 467 END DO ! jl 468 469 DO jl = 1, jpl 335 v_s (ji,jj,jl) = MAX( 0._wp, v_s (ji,jj,jl) ) 336 v_i (ji,jj,jl) = MAX( 0._wp, v_i (ji,jj,jl) ) 337 smv_i(ji,jj,jl) = MAX( 0._wp, smv_i(ji,jj,jl) ) 338 oa_i (ji,jj,jl) = MAX( 0._wp, oa_i (ji,jj,jl) ) 339 a_i (ji,jj,jl) = MAX( 0._wp, a_i (ji,jj,jl) ) 340 e_s (ji,jj,1,jl) = MAX( 0._wp, e_s (ji,jj,1,jl) ) 341 END DO 342 END DO 343 470 344 DO jk = 1, nlay_i 471 345 DO jj = 1, jpj 472 346 DO ji = 1, jpi 473 ! Ice heat content 474 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 475 ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval ) 476 e_i(ji,jj,jk,jl) = zindic * ze 477 END DO !ji 478 END DO ! jj 479 END DO ! jk 480 END DO ! jl 481 482 483 ! --- agglomerate variables (clem) ----------------- 484 vt_i (:,:) = 0._wp 485 vt_s (:,:) = 0._wp 486 at_i (:,:) = 0._wp 487 ! 488 DO jl = 1, jpl 347 e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) ) 348 END DO 349 END DO 350 END DO 351 END DO 352 !!gm & cr 353 354 ! --- diags --- 489 355 DO jj = 1, jpj 490 356 DO ji = 1, jpi 491 ! 492 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 493 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 494 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 495 ! 496 zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) ) 497 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda ! ice thickness 498 END DO 499 END DO 500 END DO 357 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 358 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 359 360 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 361 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 362 diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 363 END DO 364 END DO 365 366 ! zap small areas 367 CALL lim_var_zapsmall 368 369 !--- Thickness correction in case too high -------------------------------------------------------- 370 DO jl = 1, jpl 371 DO jj = 1, jpj 372 DO ji = 1, jpi 373 374 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 375 376 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 377 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 378 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 379 380 zvi = v_i (ji,jj,jl) 381 zvs = v_s (ji,jj,jl) 382 zsmv = smv_i(ji,jj,jl) 383 zes = e_s (ji,jj,1,jl) 384 zei = SUM( e_i(ji,jj,1:nlay_i,jl) ) 385 386 zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 387 388 IF ( ( zdv > 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 389 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 390 391 rswitch = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 392 a_i(ji,jj,jl) = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 ) 393 394 ! small correction due to *rswitch for a_i 395 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl) 396 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl) 397 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 398 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 399 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 400 401 ! Update mass fluxes 402 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 403 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 404 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 405 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 406 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 407 408 ENDIF 409 410 ENDIF 411 412 END DO 413 END DO 414 END DO 415 ! ------------------------------------------------- 416 417 !-------------------------------------- 418 ! Impose a_i < amax in mono-category 419 !-------------------------------------- 420 ! 421 IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2 422 DO jj = 1, jpj 423 DO ji = 1, jpi 424 a_i(ji,jj,1) = MIN( a_i(ji,jj,1), rn_amax ) 425 END DO 426 END DO 427 ENDIF 428 429 ! --- agglomerate variables ----------------- 430 vt_i (:,:) = 0._wp 431 vt_s (:,:) = 0._wp 432 at_i (:,:) = 0._wp 433 DO jl = 1, jpl 434 DO jj = 1, jpj 435 DO ji = 1, jpi 436 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) 437 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) 438 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) 439 END DO 440 END DO 441 END DO 442 443 ! --- open water = 1 if at_i=0 -------------------------------- 444 DO jj = 1, jpj 445 DO ji = 1, jpi 446 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 447 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 448 END DO 449 END DO 450 451 ! conservation test 452 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 453 454 ENDIF 455 501 456 ! ------------------------------------------------- 502 503 504 505 ENDIF 506 507 IF(ln_ctl) THEN ! Control print 508 CALL prt_ctl_info(' ') 509 CALL prt_ctl_info(' - Cell values : ') 510 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 511 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp : cell area :') 512 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp : at_i :') 513 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp : vt_i :') 514 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp : vt_s :') 515 DO jl = 1, jpl 516 CALL prt_ctl_info(' ') 517 CALL prt_ctl_info(' - Category : ', ivar1=jl) 518 CALL prt_ctl_info(' ~~~~~~~~~~') 519 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_trp : a_i : ') 520 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_trp : ht_i : ') 521 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_trp : ht_s : ') 522 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_trp : v_i : ') 523 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_trp : v_s : ') 524 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_trp : e_s : ') 525 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_trp : t_su : ') 526 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_trp : t_snow : ') 527 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_trp : sm_i : ') 528 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_trp : smv_i : ') 529 DO jk = 1, nlay_i 530 CALL prt_ctl_info(' ') 531 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 532 CALL prt_ctl_info(' ~~~~~~~') 533 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp : t_i : ') 534 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp : e_i : ') 535 END DO 536 END DO 537 ENDIF 538 ! ------------------------------- 539 !- check conservation (C Rousset) 540 IF( ln_limdiahsb ) THEN 541 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 542 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 543 544 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 545 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 546 547 zchk_vmin = glob_min(v_i) 548 zchk_amax = glob_max(SUM(a_i,dim=3)) 549 zchk_amin = glob_min(a_i) 550 zchk_umax = glob_max(SQRT(u_ice**2 + v_ice**2)) 551 552 IF(lwp) THEN 553 IF ( ABS( zchk_v_i ) > 1.e-5 ) THEN 554 WRITE(numout,*) 'violation volume [m3/day] (limtrp) = ',(zchk_v_i * rday) 555 WRITE(numout,*) 'u_ice max [m/s] (limtrp) = ',zchk_umax 556 WRITE(numout,*) 'number of time steps (limtrp) =',kt 557 ENDIF 558 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday) 559 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limtrp) = ',(zchk_vmin * 1.e-3) 560 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limtrp) = ',zchk_amin 561 ENDIF 562 ENDIF 563 !- check conservation (C Rousset) 564 ! ------------------------------- 457 ! control prints 458 ! ------------------------------------------------- 459 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 565 460 ! 566 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow)567 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )568 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e)569 570 CALL wrk_dealloc( jpi,jpj,jpl, zaiold, zhimax ) ! clem461 CALL wrk_dealloc( jpi,jpj, zsm, zatold, zeiold, zesold ) 462 CALL wrk_dealloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 463 CALL wrk_dealloc( jpi,jpj,1, z0opw ) 464 CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 465 CALL wrk_dealloc( jpi,jpj,jpl, zviold, zvsold, zhimax, zsmvold ) 571 466 ! 572 467 IF( nn_timing == 1 ) CALL timing_stop('limtrp') 468 573 469 END SUBROUTINE lim_trp 574 470 … … 581 477 END SUBROUTINE lim_trp 582 478 #endif 583 584 479 !!====================================================================== 585 480 END MODULE limtrp
Note: See TracChangeset
for help on using the changeset viewer.
