Changeset 5313 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
- Timestamp:
- 2015-05-29T11:46:03+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r5312 r5313 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(:,:) :: z ui_u, zvi_v, zsm, zs0at, zs0ow67 REAL(wp), POINTER, DIMENSION(:,:,:) :: z s0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi68 REAL(wp), POINTER, DIMENSION(:,:,: ,:) :: zs0e69 REAL(wp), POINTER, DIMENSION(:,:,: ) :: zviold, zvsold ! old ice volume...70 REAL(wp), POINTER, DIMENSION(:,:,:) :: z aiold, zhimax ! old ice concentration and thickness71 REAL(wp), POINTER, DIMENSION(:,: ) :: zeiold, zesold ! old enthalpies72 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei73 !74 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b70 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 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, 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 ) 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) ------------------------------- 107 CALL lim_var_glo2eqv 108 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 109 124 !--------------------------------------------------------------------- 110 ! Record max of the surrounding ice thicknesses for correction in limupdate125 ! Record max of the surrounding ice thicknesses for correction 111 126 ! in case advection creates ice too thick. 112 127 !--------------------------------------------------------------------- 113 zhimax(:,:,:) = ht_i(:,:,:) 128 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 114 129 DO jl = 1, jpl 115 130 DO jj = 2, jpjm1 116 131 DO ji = 2, jpim1 117 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) ) 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) ) 122 133 END DO 123 134 END DO … … 125 136 END DO 126 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 127 161 !------------------------- 128 162 ! transported fields 129 163 !------------------------- 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' 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 160 177 161 178 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, z s0ow (:,:), sxopw(:,:), &166 & 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(:,:) ) 167 184 DO jl = 1, jpl 168 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 --- 169 186 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 170 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), & 171 188 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 172 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 --- 173 190 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 174 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), & 175 192 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 176 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 --- 177 194 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 178 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), & 179 196 & 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 ---197 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 181 198 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 182 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), & 183 200 & 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 ---201 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 185 202 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 186 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), & 187 204 & 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 ---205 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 189 206 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 190 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), & 191 208 & 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 , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,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), & 194 211 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 195 212 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0e(:,:,jk,jl), sxe (:,:,jk,jl), &213 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 197 214 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 198 215 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) … … 201 218 END DO 202 219 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, z s0ow (:,:), sxopw(:,:), &207 & 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(:,:) ) 208 225 DO jl = 1, jpl 209 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 --- 210 227 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 211 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), & 212 229 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 213 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 --- 214 231 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 215 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), & 216 233 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 217 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 --- 218 235 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 219 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), & 220 237 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 221 222 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 --- 223 239 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 224 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), & 225 241 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 226 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 --- 227 243 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 228 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), & 229 245 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 230 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 --- 231 247 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 232 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), & 233 249 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 234 250 DO jk = 1, nlay_i !--- ice heat contents --- 235 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &251 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 236 252 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 237 253 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0e(:,:,jk,jl), sxe (:,:,jk,jl), &254 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 239 255 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 240 256 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) … … 247 263 ! Recover the properties from their contents 248 264 !------------------------------------------- 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 ! 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(:,:) 273 DO jk = 1, nlay_i 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) 257 281 END DO 258 282 259 283 !------------------------------------------------------------------------------! 260 ! 4)Diffusion of Ice fields284 ! Diffusion of Ice fields 261 285 !------------------------------------------------------------------------------! 262 286 287 ! 263 288 !-------------------------------- 264 289 ! diffusion of open water area 265 290 !-------------------------------- 266 zs0at(:,:) = zs0a(:,:,1) ! total ice fraction267 DO jl = 2, jpl268 zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)269 END DO270 !271 291 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points 272 292 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 273 293 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)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) 278 298 END DO 279 299 END DO 280 300 ! 281 CALL lim_hdf( zs0ow (:,:) ) ! Diffusion301 CALL lim_hdf( ato_i (:,:) ) 282 302 283 303 !------------------------------------ … … 288 308 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 289 309 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) )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) ) 303 323 DO jk = 1, nlay_i 304 CALL lim_hdf( zs0e(:,:,jk,jl) )324 CALL lim_hdf( e_i(:,:,jk,jl) ) 305 325 END DO 306 326 END DO 307 327 308 328 !------------------------------------------------------------------------------! 309 ! 5) Update andlimit ice properties after transport329 ! limit ice properties after transport 310 330 !------------------------------------------------------------------------------! 311 312 !-------------------------------------------------- 313 ! 5.1) Recover mean values over the grid squares. 314 !-------------------------------------------------- 315 zs0at(:,:) = 0._wp 331 !!gm & cr : MAX should not be active if adv scheme is positive ! 316 332 DO jl = 1, jpl 317 333 DO jj = 1, jpj 318 334 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 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 364 344 DO jk = 1, nlay_i 365 345 DO jj = 1, jpj 366 346 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 fluxes 371 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 372 END DO !ji 373 END DO ! jj 374 END DO ! jk 375 END DO ! jl 376 377 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 378 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 -------------------------------------------------------- 379 370 DO jl = 1, jpl 380 371 DO jj = 1, jpj … … 382 373 383 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 384 380 zvi = v_i (ji,jj,jl) 385 381 zvs = v_s (ji,jj,jl) … … 387 383 zes = e_s (ji,jj,1,jl) 388 384 zei = SUM( e_i(ji,jj,1:nlay_i,jl) ) 389 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 390 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 391 392 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. & 394 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 395 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 396 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 397 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 398 ELSE 399 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 400 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 401 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 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 402 408 ENDIF 403 409 404 ! small correction due to *rswitch for a_i405 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl)406 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl)407 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)408 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)409 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)410 411 ! Update mass fluxes412 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice413 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice414 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice415 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 <0417 410 ENDIF 411 418 412 END DO 419 413 END DO 420 414 END DO 421 415 ! ------------------------------------------------- 422 423 ! --- diags --- 424 DO jj = 1, jpj 425 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 431 END DO 432 END DO 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 433 428 434 429 ! --- agglomerate variables ----------------- … … 436 431 vt_s (:,:) = 0._wp 437 432 at_i (:,:) = 0._wp 438 !439 433 DO jl = 1, jpl 440 434 DO jj = 1, jpj 441 435 DO ji = 1, jpi 442 ! 443 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 444 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 445 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 446 END DO 447 END DO 448 END DO 449 ! ------------------------------------------------- 450 451 ! 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 -------------------------------- 452 444 DO jj = 1, jpj 453 445 DO ji = 1, jpi 454 ! open water = 1 if at_i=0455 446 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 456 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj)447 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 457 448 END DO 458 449 END DO … … 463 454 ENDIF 464 455 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 456 ! ------------------------------------------------- 457 ! control prints 458 ! ------------------------------------------------- 459 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 496 460 ! 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 ) ! 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 ) 502 466 ! 503 467 IF( nn_timing == 1 ) CALL timing_stop('limtrp') 468 504 469 END SUBROUTINE lim_trp 505 470 … … 512 477 END SUBROUTINE lim_trp 513 478 #endif 514 515 479 !!====================================================================== 516 480 END MODULE limtrp
Note: See TracChangeset
for help on using the changeset viewer.