Changeset 924 for codes/icosagcm/devel/src/dissip/dissip_gcm.f90
- Timestamp:
- 06/19/19 17:22:32 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dissip/dissip_gcm.f90
r856 r924 26 26 REAL,ALLOCATABLE,SAVE :: tau_divgrad(:) 27 27 !$OMP THREADPRIVATE(tau_divgrad) 28 28 29 29 REAL,SAVE :: cgraddiv 30 30 !$OMP THREADPRIVATE(cgraddiv) … … 33 33 REAL,SAVE :: cdivgrad 34 34 !$OMP THREADPRIVATE(cdivgrad) 35 35 36 36 INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear 37 37 !$OMP THREADPRIVATE(rayleigh_friction_type) … … 126 126 CALL dissip_constants 127 127 CALL dissip_profile 128 CALL dissip_timescale 128 129 END IF 129 130 … … 236 237 SUBROUTINE dissip_profile 237 238 USE disvert_mod 238 INTEGER :: l 239 REAL(rstd) :: mintau, zz, zvert, fact 240 fact=2 239 ! parameters used by the various profiles 240 ! IF planet_type == earth 241 ! IF dissip_vert_prof == 0 242 ! none 243 ! IF dissip_vert_prof == 1 244 ! dissip_zref, dissip_deltaz, dissip_factz 245 ! IF planet_type == other 246 ! IF dissip_vert_prof == 0 247 ! dissip_fac_mid 248 ! + dissip_deltaz, dissip_hdelta, dissip_fac_up, dissip_pupstart IF ok_strato 249 ! IF dissip_vert_prof == 1 250 ! fac_mid, fac_up, startalt, delta => middle (hardcoded), scaleheight 251 ! 252 253 REAL(rstd), PARAMETER :: fact=2., & 254 fac_mid=3., & ! coefficient for lower/middle atmosphere 255 fac_up=30., & ! coefficient for upper atmosphere 256 startalt=70., & ! altitude (in km) for the transition from middle to upper atm. 257 delta=30., & ! Size (in km) of the transition region between middle/upper atm. 258 middle=startalt+delta/2. 259 REAL(rstd) :: dissip_zref, dissip_deltaz, dissip_factz, & 260 dissip_hdelta, dissip_fac_up, dissip_fac_mid, dissip_pupstart, & 261 scaleheight, & 262 zz, pseudoz, pup, & 263 sigma(llm), zvert(llm) 264 CHARACTER(LEN=255) :: planet_type ! earth, other, other_strato ; other_strato = other + ok_strato 265 LOGICAL :: ok_strato 266 INTEGER :: l, dissip_vert_prof 267 268 ! select vertical profile of horizontal dissipation coefficients, see also [781] 269 270 planet_type='earth' 271 CALL getin('dissip_planet_type', planet_type) 272 SELECT CASE(planet_type) 273 CASE('earth','other') 274 ok_strato=.FALSE. 275 CASE('other_strato') 276 planet_type='other' 277 ok_strato=.TRUE. 278 CASE DEFAULT 279 STOP "Invalid value of dissip_planet_type, valid values are <earth>, <other>, <other_strato>" 280 END SELECT 281 282 dissip_vert_prof = 0 283 CALL getin('dissip_vert_prof',dissip_vert_prof) 284 285 SELECT CASE(dissip_vert_prof) 286 CASE(0) 287 IF(TRIM(planet_type)=='other') THEN 288 ! Default values given below are for a Venus-like planet 289 CALL getin('dissip_fac_mid',dissip_fac_mid ) 290 dissip_fac_mid=2. 291 IF(ok_strato) THEN 292 dissip_fac_up=50. 293 ! deltaz et hdelta in km 294 dissip_deltaz=30. 295 dissip_hdelta=5. 296 ! pupstart in Pa 297 dissip_pupstart=1.e4 298 CALL getin('dissip_deltaz',dissip_deltaz ) 299 CALL getin('dissip_hdelta',dissip_hdelta ) 300 CALL getin('dissip_fac_up',dissip_fac_up ) 301 CALL getin('dissip_pupstart',dissip_pupstart) 302 END IF 303 END IF 304 CASE(1) 305 IF(TRIM(planet_type)=='earth') THEN 306 dissip_zref = 30. 307 dissip_deltaz = 10. 308 dissip_factz = 4. 309 CALL getin('dissip_zref',dissip_zref ) 310 CALL getin('dissip_deltaz',dissip_deltaz ) 311 CALL getin('dissip_factz',dissip_factz ) 312 ELSE 313 ! fac_mid, fac_up, startalt, delta => middle are hardcoded 314 CALL getin('dissip_scaleheight', scaleheight) 315 END IF 316 CASE DEFAULT 317 STOP 'Invalid value of dissip_vert_prof : valid values are 0,1' 318 END SELECT 319 320 IF(ap_bp_present) THEN 321 sigma(:) = preff/presnivs(:) 322 ELSE 323 sigma(:) = 1. 324 END IF 325 326 SELECT CASE(TRIM(planet_type)) 327 CASE('earth') 328 DO l=1,llm 329 IF(dissip_vert_prof == 1) THEN 330 pseudoz=8.*LOG(sigma(l)) 331 zvert(l)=1+ & 332 (TANH((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. & 333 *(dissip_factz-1.) 334 ELSE 335 zz = 1.-sigma(l) 336 zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) 337 END IF 338 END DO 339 CASE('other') 340 SELECT CASE(dissip_vert_prof) 341 CASE(0) 342 ! First step: going from 1 to dissip_fac_mid (in gcm.def) 343 !============ 344 DO l=1,llm 345 zz = 1. - sigma(l) 346 zvert(l) = dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 347 END DO 348 ! 349 ! Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 350 !========================== 351 ! Utilisation de la fonction tangente hyperbolique pour augmenter 352 ! arbitrairement la dissipation et donc la stabilite du modele a 353 ! partir d'une certaine altitude. 354 ! 355 ! Le facteur multiplicatif de basse atmosphere etant deja pris 356 ! en compte, il faut diviser le facteur multiplicatif de haute 357 ! atmosphere par celui-ci. 358 359 IF(ok_strato) THEN 360 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 361 DO l=1,llm 362 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) & 363 *(1.-(0.5*(1+tanh(-6./dissip_deltaz* & 364 (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 365 END DO 366 END IF 367 CASE(1) 368 DO l=1,llm 369 zz = 1. - sigma(l) 370 zvert(l)= fac_mid -( fac_mid-1.)/( 1.+zz*zz ) 371 zvert(l)= zvert(l)*(1.0+((fac_up/fac_mid-1)* & 372 (1.-(0.5*(1+tanh(-6./ & 373 delta*(scaleheight*(-log(presnivs(l)/preff))-middle))))) & 374 )) 375 END DO 376 END SELECT 377 END SELECT 378 241 379 DO l=1,llm 242 IF(ap_bp_present) THEN ! height-dependent dissipation 243 zz= 1. - preff/presnivs(l) 244 ELSE 245 zz = 0. 246 END IF 247 zvert=fact-(fact-1)/(1+zz*zz) 248 tau_graddiv(l) = tau_graddiv(l)/zvert 249 tau_gradrot(l) = tau_gradrot(l)/zvert 250 tau_divgrad(l) = tau_divgrad(l)/zvert 380 tau_graddiv(l) = tau_graddiv(l)/zvert(l) 381 tau_gradrot(l) = tau_gradrot(l)/zvert(l) 382 tau_divgrad(l) = tau_divgrad(l)/zvert(l) 251 383 END DO 252 384 END SUBROUTINE dissip_profile 385 386 SUBROUTINE dissip_timescale 387 INTEGER :: l 388 REAL(rstd) :: mintau 253 389 mintau=tau_graddiv(1) 254 390 DO l=1,llm … … 276 412 ENDIF 277 413 278 END SUBROUTINE dissip_ profile279 414 END SUBROUTINE dissip_timescale 415 280 416 SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due) 281 417 TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
Note: See TracChangeset
for help on using the changeset viewer.