- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4688 r6225 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 39 REAL(wp) :: epsi10 = 1.e-10_wp 40 REAL(wp) :: epsi20 = 1.e-20_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) :: zcfl , zusnit ! - - 70 REAL(wp) :: zsal , zage ! - - 67 REAL(wp) :: zcfl , zusnit ! - - 68 CHARACTER(len=80) :: cltmp 71 69 ! 72 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 73 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 74 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 75 ! mass and salt flux (clem) 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 78 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 79 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 80 ! 81 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 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 82 79 !!--------------------------------------------------------------------- 83 80 IF( nn_timing == 1 ) CALL timing_start('limtrp') 84 81 85 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )86 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )87 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e)88 89 CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem82 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 ) 90 87 91 88 IF( numit == nstart .AND. lwp ) THEN … … 95 92 ENDIF 96 93 WRITE(numout,*) '~~~~~~~~~~~~' 94 ncfl = 0 ! nb of time step with CFL > 1/2 97 95 ENDIF 96 97 zsm(:,:) = e1e2t(:,:) 98 98 99 zsm(:,:) = area(:,:)100 101 99 ! !-------------------------------------! 102 100 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! … … 104 102 105 103 ! conservation test 106 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)107 108 ! 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 109 107 zviold(:,:,:) = v_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. (clem) ------------------------------- 114 CALL lim_var_glo2eqv 115 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 116 124 !--------------------------------------------------------------------- 117 ! Record max of the surrounding ice thicknesses for correction in limupdate125 ! Record max of the surrounding ice thicknesses for correction 118 126 ! in case advection creates ice too thick. 119 127 !--------------------------------------------------------------------- 120 zhimax(:,:,:) = ht_i(:,:,:) 128 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 121 129 DO jl = 1, jpl 122 130 DO jj = 2, jpjm1 123 131 DO ji = 2, jpim1 124 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) ) 125 !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) & 126 ! & + ht_i(ji-1,jj ,jl) * tmask(ji-1,jj ,1) + ht_i(ji ,jj-1,jl) * tmask(ji ,jj-1,1) & 127 ! & + ht_i(ji+1,jj ,jl) * tmask(ji+1,jj ,1) + ht_i(ji ,jj+1,jl) * tmask(ji ,jj+1,1) & 128 ! & + 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) ) 129 133 END DO 130 134 END DO … … 132 136 END DO 133 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 134 161 !------------------------- 135 162 ! transported fields 136 163 !------------------------- 137 ! Snow vol, ice vol, salt and age contents, area 138 zs0ow(:,:) = ato_i(:,:) * area(:,:) ! Open water area 139 DO jl = 1, jpl 140 zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume 141 zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume 142 zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area 143 zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content 144 zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content 145 zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content 146 zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content 147 END DO 148 149 !-------------------------- 150 ! Advection of Ice fields (Prather scheme) 151 !-------------------------- 152 ! If ice drift field is too fast, use an appropriate time step for advection. 153 ! CFL test for stability 154 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) ) 155 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 156 IF(lk_mpp ) CALL mpp_max( zcfl ) 157 !!gm more readability: 158 ! IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 159 ! ELSE ; initad = 1 ; zusnit = 1.0_wp 160 ! ENDIF 161 !!gm end 162 initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) ) 163 zusnit = 1.0 / REAL( initad ) 164 IF( zcfl > 0.5 .AND. lwp ) & 165 WRITE(numout,*) 'lim_trp : CFL violation at day ', nday, ', cfl = ', zcfl, & 166 & ': the ice time stepping is split in two' 164 z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:) ! Open water area 165 DO jl = 1, jpl 166 z0snw (:,:,jl) = v_s (:,:, jl) * e1e2t(:,:) ! Snow volume 167 z0ice(:,:,jl) = v_i (:,:, jl) * e1e2t(:,:) ! Ice volume 168 z0ai (:,:,jl) = a_i (:,:, jl) * e1e2t(:,:) ! Ice area 169 z0smi (:,:,jl) = smv_i(:,:, jl) * e1e2t(:,:) ! Salt content 170 z0oi (:,:,jl) = oa_i (:,:, jl) * e1e2t(:,:) ! Age content 171 z0es (:,:,jl) = e_s (:,:,1,jl) * e1e2t(:,:) ! Snow heat content 172 DO jk = 1, nlay_i 173 z0ei (:,:,jk,jl) = e_i (:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 174 END DO 175 END DO 176 167 177 168 178 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 169 DO j k = 1,initad170 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area171 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )172 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0ow (:,:), sxopw(:,:), &173 & 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(:,:) ) 174 184 DO jl = 1, jpl 175 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---185 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 176 186 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 177 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0ice(:,:,jl), sxice(:,:,jl), &187 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & 178 188 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 179 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn(:,:,jl), sxsn (:,:,jl), & !--- snow volume ---189 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 180 190 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 181 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0sn(:,:,jl), sxsn (:,:,jl), &191 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & 182 192 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 183 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm(:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---193 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 184 194 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 185 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0sm(:,:,jl), sxsal(:,:,jl), &195 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & 186 196 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 187 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), &!--- ice age ---197 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 188 198 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 189 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0oi(:,:,jl), sxage(:,:,jl), &199 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & 190 200 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 191 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), &!--- ice concentrations ---201 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 192 202 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 193 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0a(:,:,jl), sxa (:,:,jl), &203 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & 194 204 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 195 CALL lim_adv_x( zusnit, u_ice, 1._wp , 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 --- 196 206 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 197 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0c0(:,:,jl), sxc0 (:,:,jl), &207 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & 198 208 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 199 DO layer = 1, nlay_i!--- ice heat contents ---200 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &201 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &202 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )203 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0e(:,:,layer,jl), sxe (:,:,layer,jl), &204 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &205 & 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) ) 206 216 END DO 207 217 END DO 208 218 END DO 209 219 ELSE 210 DO j k= 1, initad211 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area212 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )213 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0ow (:,:), sxopw(:,:), &214 & 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(:,:) ) 215 225 DO jl = 1, jpl 216 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---226 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 217 227 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 218 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0ice(:,:,jl), sxice(:,:,jl), &228 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & 219 229 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 220 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn(:,:,jl), sxsn (:,:,jl), & !--- snow volume ---230 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 221 231 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 222 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0sn(:,:,jl), sxsn (:,:,jl), &232 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & 223 233 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 224 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm(:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---234 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 225 235 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 226 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0sm(:,:,jl), sxsal(:,:,jl), &236 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & 227 237 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 228 229 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 238 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 230 239 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 231 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0oi(:,:,jl), sxage(:,:,jl), &240 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & 232 241 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 233 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a(:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---242 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 234 243 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 235 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0a(:,:,jl), sxa (:,:,jl), &244 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & 236 245 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 237 CALL lim_adv_y( zusnit, v_ice, 1._wp , 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 --- 238 247 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 239 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0c0(:,:,jl), sxc0 (:,:,jl), &248 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & 240 249 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 241 DO layer= 1, nlay_i !--- ice heat contents ---242 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &243 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &244 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )245 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0e(:,:,layer,jl), sxe (:,:,layer,jl), &246 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &247 & 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) ) 248 257 END DO 249 258 END DO … … 254 263 ! Recover the properties from their contents 255 264 !------------------------------------------- 256 zs0ow(:,:) = zs0ow(:,:) / area(:,:) 257 DO jl = 1, jpl 258 zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:) 259 zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:) 260 zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:) 261 zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 262 zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:) 263 ! 265 ato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) 266 DO jl = 1, jpl 267 v_i (:,:, jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) 268 v_s (:,:, jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) 269 smv_i(:,:, jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) 270 oa_i (:,:, jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) 271 a_i (:,:, jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) 272 e_s (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:) 273 DO jk = 1, nlay_i 274 e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) 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) 264 281 END DO 265 282 266 283 !------------------------------------------------------------------------------! 267 ! 4)Diffusion of Ice fields284 ! Diffusion of Ice fields 268 285 !------------------------------------------------------------------------------! 269 286 287 ! 270 288 !-------------------------------- 271 289 ! diffusion of open water area 272 290 !-------------------------------- 273 zs0at(:,:) = zs0a(:,:,1) ! total ice fraction274 DO jl = 2, jpl275 zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)276 END DO277 !278 291 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points 279 292 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 280 293 DO ji = 1 , fs_jpim1 ! vector opt. 281 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0at(ji ,jj) ) ) ) &282 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)283 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0at(ji,jj ) ) ) ) &284 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- 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) 285 298 END DO 286 299 END DO 287 300 ! 288 CALL lim_hdf( zs0ow (:,:) ) ! Diffusion301 CALL lim_hdf( ato_i (:,:) ) 289 302 290 303 !------------------------------------ … … 295 308 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 296 309 DO ji = 1 , fs_jpim1 ! vector opt. 297 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0a(ji ,jj,jl) ) ) ) &298 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)299 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - zs0a(ji,jj ,jl) ) ) ) &300 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)301 END DO 302 END DO 303 304 CALL lim_hdf( zs0ice (:,:,jl) )305 CALL lim_hdf( zs0sn (:,:,jl) )306 CALL lim_hdf( zs0sm (:,:,jl) )307 CALL lim_hdf( zs0oi (:,:,jl) )308 CALL lim_hdf( zs0a (:,:,jl) )309 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) ) 310 323 DO jk = 1, nlay_i 311 CALL lim_hdf( zs0e(:,:,jk,jl) )324 CALL lim_hdf( e_i(:,:,jk,jl) ) 312 325 END DO 313 326 END DO 314 327 315 328 !------------------------------------------------------------------------------! 316 ! 5) Update andlimit ice properties after transport329 ! limit ice properties after transport 317 330 !------------------------------------------------------------------------------! 318 319 !-------------------------------------------------- 320 ! 5.1) Recover mean values over the grid squares. 321 !-------------------------------------------------- 322 zs0at(:,:) = 0._wp 331 !!gm & cr : MAX should not be active if adv scheme is positive ! 323 332 DO jl = 1, jpl 324 333 DO jj = 1, jpj 325 334 DO ji = 1, jpi 326 zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) ) 327 zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) ) 328 zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) ) 329 zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) ) 330 zs0a (ji,jj,jl) = MAX( 0._wp, zs0a (ji,jj,jl) ) 331 zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) ) 332 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 333 END DO 334 END DO 335 END DO 336 337 !--------------------------------------------------------- 338 ! 5.2) Update and mask variables 339 !--------------------------------------------------------- 340 DO jl = 1, jpl 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 344 345 zvi = zs0ice(ji,jj,jl) 346 zvs = zs0sn (ji,jj,jl) 347 zes = zs0c0 (ji,jj,jl) 348 zsmv = zs0sm (ji,jj,jl) 349 ! 350 ! Remove very small areas 351 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 352 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 353 a_i(ji,jj,jl) = zindb * zs0a (ji,jj,jl) 354 e_s(ji,jj,1,jl) = zindb * zs0c0 (ji,jj,jl) 355 ! Ice salinity and age 356 IF( num_sal == 2 ) THEN 357 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 358 ENDIF 359 oa_i(ji,jj,jl) = MAX( zindb * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 360 361 ! Update fluxes 362 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 363 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 364 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 365 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 366 END DO 367 END DO 368 END DO 369 370 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 371 344 DO jk = 1, nlay_i 372 345 DO jj = 1, jpj 373 346 DO ji = 1, jpi 374 zindb = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 375 zei = zs0e(ji,jj,jk,jl) 376 e_i(ji,jj,jk,jl) = zindb * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 377 ! Update fluxes 378 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 <0 379 END DO !ji 380 END DO ! jj 381 END DO ! jk 382 END DO ! jl 383 384 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 385 CALL lim_var_glo2eqv 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 --- 355 DO jj = 1, jpj 356 DO ji = 1, jpi 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 -------------------------------------------------------- 386 370 DO jl = 1, jpl 387 371 DO jj = 1, jpj … … 389 373 390 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 391 380 zvi = v_i (ji,jj,jl) 392 381 zvs = v_s (ji,jj,jl) 393 382 zsmv = smv_i(ji,jj,jl) 394 383 zes = e_s (ji,jj,1,jl) 395 zei = SUM( e_i(ji,jj,:,jl) ) 396 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 397 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 398 399 zindh = 1._wp 400 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 401 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 402 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 403 zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 404 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 405 ELSE 406 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 407 zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 408 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 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 409 408 ENDIF 410 409 411 ! small correction due to *zindh for a_i412 v_i (ji,jj,jl) = zindh * v_i (ji,jj,jl)413 v_s (ji,jj,jl) = zindh * v_s (ji,jj,jl)414 smv_i(ji,jj,jl) = zindh * smv_i(ji,jj,jl)415 e_s(ji,jj,1,jl) = zindh * e_s(ji,jj,1,jl)416 e_i(ji,jj,:,jl) = zindh * e_i(ji,jj,:,jl)417 418 ! Update mass fluxes419 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice420 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice421 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice422 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 <0423 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,:,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0424 425 410 ENDIF 426 411 427 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice428 diag_trp_vs(ji,jj) = diag_trp_vs(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * r1_rdtice429 430 412 END DO 431 413 END DO 432 414 END DO 433 415 ! ------------------------------------------------- 434 435 ! --- diags --- 436 DO jj = 1, jpj 437 DO ji = 1, jpi 438 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 439 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 440 END DO 441 END DO 442 443 ! --- agglomerate variables (clem) ----------------- 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 ----------------- 444 430 vt_i (:,:) = 0._wp 445 431 vt_s (:,:) = 0._wp 446 432 at_i (:,:) = 0._wp 447 !448 433 DO jl = 1, jpl 449 434 DO jj = 1, jpj 450 435 DO ji = 1, jpi 451 ! 452 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 453 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 454 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 455 END DO 456 END DO 457 END DO 458 ! ------------------------------------------------- 459 460 ! open water 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 -------------------------------- 461 444 DO jj = 1, jpj 462 445 DO ji = 1, jpi 463 ! open water = 1 if at_i=0 464 zindb = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 465 ato_i(ji,jj) = zindb + (1._wp - zindb ) * zs0ow(ji,jj) 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) 466 448 END DO 467 449 END DO … … 472 454 ENDIF 473 455 474 IF(ln_ctl) THEN ! Control print 475 CALL prt_ctl_info(' ') 476 CALL prt_ctl_info(' - Cell values : ') 477 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 478 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp : cell area :') 479 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp : at_i :') 480 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp : vt_i :') 481 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp : vt_s :') 482 DO jl = 1, jpl 483 CALL prt_ctl_info(' ') 484 CALL prt_ctl_info(' - Category : ', ivar1=jl) 485 CALL prt_ctl_info(' ~~~~~~~~~~') 486 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_trp : a_i : ') 487 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_trp : ht_i : ') 488 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_trp : ht_s : ') 489 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_trp : v_i : ') 490 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_trp : v_s : ') 491 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_trp : e_s : ') 492 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_trp : t_su : ') 493 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_trp : t_snow : ') 494 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_trp : sm_i : ') 495 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_trp : smv_i : ') 496 DO jk = 1, nlay_i 497 CALL prt_ctl_info(' ') 498 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 499 CALL prt_ctl_info(' ~~~~~~~') 500 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp : t_i : ') 501 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp : e_i : ') 502 END DO 503 END DO 504 ENDIF 456 ! ------------------------------------------------- 457 ! control prints 458 ! ------------------------------------------------- 459 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 505 460 ! 506 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )507 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )508 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e)509 510 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, 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 ) 511 466 ! 512 467 IF( nn_timing == 1 ) CALL timing_stop('limtrp') 468 513 469 END SUBROUTINE lim_trp 514 470 … … 521 477 END SUBROUTINE lim_trp 522 478 #endif 523 524 479 !!====================================================================== 525 480 END MODULE limtrp
Note: See TracChangeset
for help on using the changeset viewer.
