- Timestamp:
- 2017-09-12T20:46:13+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
r8516 r8517 37 37 PUBLIC ice_adv_init ! called by icestp 38 38 39 INTEGER :: nice_adv ! choice of the type of advection scheme 40 ! ! associated indices: 41 INTEGER, PARAMETER :: np_advPRA = 1 ! Prather scheme 42 INTEGER, PARAMETER :: np_advUMx = 2 ! Ultimate-Macho scheme 43 39 44 !! * Substitution 40 45 # include "vectopt_loop_substitute.h90" … … 61 66 !!---------------------------------------------------------------------- 62 67 INTEGER, INTENT(in) :: kt ! number of iteration 63 !64 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices65 INTEGER :: initad ! number of sub-timestep for the advection66 REAL(wp) :: zcfl , zusnit ! - -67 CHARACTER(len=80) :: cltmp68 !69 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b70 REAL(wp) :: zdv71 REAL(wp), DIMENSION(jpi,jpj) :: zatold, zeiold, zesold, zsmvold72 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhimax, zviold, zvsold73 68 !!--------------------------------------------------------------------- 74 69 IF( nn_timing == 1 ) CALL timing_start('iceadv') … … 80 75 ENDIF 81 76 82 CALL ice_var_agg( 1 ) ! integrated values + ato_i83 84 77 ! conservation test 85 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 86 87 ! store old values for diag 88 zviold (:,:,:) = v_i(:,:,:) 89 zvsold (:,:,:) = v_s(:,:,:) 90 zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 ) 91 zeiold (:,:) = et_i(:,:) 92 zesold (:,:) = et_s(:,:) 93 94 ! Thickness correction init. 95 zatold(:,:) = at_i 96 WHERE( a_i(:,:,:) >= epsi20 ) 97 ht_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 98 ht_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:) 99 ELSEWHERE 100 ht_i(:,:,:) = 0._wp 101 ht_s(:,:,:) = 0._wp 102 END WHERE 103 104 ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick 105 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 106 DO jl = 1, jpl 107 DO jj = 2, jpjm1 108 DO ji = 2, jpim1 109 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 110 zhimax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) ) 111 END DO 112 END DO 113 END DO 114 CALL lbc_lnk( zhimax(:,:,:), 'T', 1. ) 115 78 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 79 116 80 !---------- 117 81 ! Advection 118 82 !---------- 119 IF( ln_adv_UMx ) THEN !-- ULTIMATE-MACHO scheme 83 SELECT CASE( nice_adv ) 84 ! !-----------------------! 85 CASE( np_advUMx ) ! ULTIMATE-MACHO scheme ! 86 ! !-----------------------! 120 87 CALL ice_adv_umx( kt, u_ice, v_ice, & 121 88 & ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 122 123 ELSEIF( ln_adv_Pra ) THEN !-- PRATHER scheme 89 ! !-----------------------! 90 CASE( np_advPRA ) ! PRATHER scheme ! 91 ! !-----------------------! 124 92 CALL ice_adv_prather( kt, u_ice, v_ice, & 125 93 & ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 126 127 ENDIF 94 END SELECT 128 95 129 ! total ice fraction130 at_i(:,:) = a_i(:,:,1)131 DO jl = 2, jpl132 at_i(:,:) = at_i(:,:) + a_i(:,:,jl)133 END DO134 135 96 !------------ 136 97 ! diagnostics 137 98 !------------ 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 diag_trp_ei (ji,jj) = ( SUM( e_i (ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 141 diag_trp_es (ji,jj) = ( SUM( e_s (ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 142 diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:) ) - zsmvold(ji,jj) ) * r1_rdtice 143 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 144 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 145 END DO 146 END DO 99 diag_trp_ei (:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice 100 diag_trp_es (:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice 101 diag_trp_smv(:,:) = SUM( smv_i(:,:,:) - smv_i_b(:,:,:) , dim=3 ) * r1_rdtice 102 diag_trp_vi (:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_rdtice 103 diag_trp_vs (:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_rdtice 147 104 IF( iom_use('icetrp') ) CALL iom_put( "icetrp" , diag_trp_vi * rday ) ! ice volume transport 148 105 IF( iom_use('snwtrp') ) CALL iom_put( "snwtrp" , diag_trp_vs * rday ) ! snw volume transport … … 150 107 IF( iom_use('deitrp') ) CALL iom_put( "deitrp" , diag_trp_ei ) ! advected ice enthalpy (W/m2) 151 108 IF( iom_use('destrp') ) CALL iom_put( "destrp" , diag_trp_es ) ! advected snw enthalpy (W/m2) 152 153 !-------------------------------------- 154 ! Thickness correction in case too high 155 !-------------------------------------- 156 IF( nn_icedyn == 2 ) THEN 157 ! 158 CALL ice_var_zapsmall !-- zap small areas 159 ! 160 DO jl = 1, jpl 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 IF ( v_i(ji,jj,jl) > 0._wp ) THEN !-- bound to zhimax 164 ! 165 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 166 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / a_i(ji,jj,jl) 167 zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 168 ! 169 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. & 170 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 171 a_i (ji,jj,jl) = ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / zhimax(ji,jj,jl) 172 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 173 ENDIF 174 ! 175 ENDIF 176 END DO 177 END DO 178 END DO 179 180 WHERE( ht_i(:,:,jpl) > hi_max(jpl) ) !-- bound ht_i to hi_max (99 m) 181 ht_i(:,:,jpl) = hi_max(jpl) 182 a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl) 183 END WHERE 184 185 IF ( nn_pnd_scheme > 0 ) THEN !-- correct pond fraction to avoid a_ip > a_i 186 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:) 187 ENDIF 188 ! 189 ENDIF 190 191 !------------------------------------------------------------ 192 ! Impose a_i < amax if no ridging/rafting or in mono-category 193 !------------------------------------------------------------ 194 IF( l_piling ) THEN !-- simple conservative piling, comparable with 1-cat models 195 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 196 DO jl = 1, jpl 197 WHERE( at_i(:,:) > epsi20 ) 198 a_i(:,:,jl) = a_i(:,:,jl) * ( 1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:) ) 199 END WHERE 200 END DO 201 ENDIF 202 203 ! agglomerate variables 204 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 205 vt_s(:,:) = SUM( v_s(:,:,:), dim=3 ) 206 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 207 208 ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i) 209 IF ( nn_pnd_scheme > 0 ) THEN 210 at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 ) 211 vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 ) 212 ENDIF 213 ! END MP 2016 214 215 ! open water = 1 if at_i=0 216 WHERE( at_i == 0._wp ) ato_i = 1._wp 217 218 ! conservation test 219 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 220 221 ! -------------- 222 ! control prints 223 ! -------------- 224 IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 225 ! 226 IF( nn_timing == 1 ) CALL timing_stop('iceadv') 109 110 ! controls 111 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 112 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ') ! prints 113 IF( nn_timing == 1 ) CALL timing_stop ('iceadv') ! timing 227 114 ! 228 115 END SUBROUTINE ice_adv … … 241 128 !! ** input : Namelist namice_adv 242 129 !!------------------------------------------------------------------- 243 INTEGER :: ios ! Local integer output status for namelist read130 INTEGER :: ios, ioptio ! Local integer output status for namelist read 244 131 !! 245 132 NAMELIST/namice_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx … … 260 147 WRITE(numout,*) '~~~~~~~~~~~~' 261 148 WRITE(numout,*) ' Namelist namice_adv' 262 WRITE(numout,*) ' advection scheme for ice transport (limtrp)' 263 WRITE(numout,*) ' type of advection scheme (Prather) ln_adv_Pra = ', ln_adv_Pra 264 WRITE(numout,*) ' type of advection scheme (Ulimate-Macho) ln_adv_UMx = ', ln_adv_UMx 265 WRITE(numout,*) ' order of the Ultimate-Macho scheme nn_UMx = ', nn_UMx 149 WRITE(numout,*) ' type of advection scheme (Prather) ln_adv_Pra = ', ln_adv_Pra 150 WRITE(numout,*) ' type of advection scheme (Ulimate-Macho) ln_adv_UMx = ', ln_adv_UMx 151 WRITE(numout,*) ' order of the Ultimate-Macho scheme nn_UMx = ', nn_UMx 266 152 ENDIF 267 153 ! 268 IF ( ( ln_adv_Pra .AND. ln_adv_UMx ) .OR. ( .NOT.ln_adv_Pra .AND. .NOT.ln_adv_UMx ) ) THEN 269 CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 270 ENDIF 154 ! !== set the choice of ice advection ==! 155 ioptio = 0 156 IF( ln_adv_Pra ) THEN ; ioptio = ioptio + 1 ; nice_adv = np_advPRA ; ENDIF 157 IF( ln_adv_UMx ) THEN ; ioptio = ioptio + 1 ; nice_adv = np_advUMx ; ENDIF 158 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 271 159 ! 272 160 END SUBROUTINE ice_adv_init
Note: See TracChangeset
for help on using the changeset viewer.
