[4449] | 1 | module atke_exchange_coeff_mod |
---|
| 2 | |
---|
| 3 | implicit none |
---|
| 4 | |
---|
| 5 | contains |
---|
| 6 | |
---|
| 7 | subroutine atke_compute_km_kh(ngrid,nlay, & |
---|
| 8 | wind_u,wind_v,temp,play,pinterf, & |
---|
| 9 | tke,Km,Kh) |
---|
| 10 | |
---|
| 11 | !======================================================================== |
---|
| 12 | ! Routine that computes turbulent Km / Kh coefficients with a |
---|
| 13 | ! 1.5 order closure scheme (TKE) with or without stationarity assumption |
---|
| 14 | ! |
---|
| 15 | ! This parameterization has been constructed in the framework of a |
---|
| 16 | ! collective and collaborative workshop, |
---|
| 17 | ! the so-called 'Atelier TKE (ATKE)' with |
---|
| 18 | ! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, |
---|
| 19 | ! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon |
---|
| 20 | ! |
---|
| 21 | ! Main assumptions of the model : |
---|
| 22 | ! (1) dry atmosphere |
---|
| 23 | ! (2) horizontal homogeneity (Dx=Dy=0.) |
---|
| 24 | !======================================================================= |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | |
---|
| 28 | USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cn, cinf, rpi |
---|
| 29 | USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, rg, rd |
---|
| 30 | |
---|
| 31 | implicit none |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | ! Declarations: |
---|
| 35 | !============= |
---|
| 36 | |
---|
| 37 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
---|
| 38 | INTEGER, INTENT(IN) :: nlay ! number of vertical index |
---|
| 39 | |
---|
| 40 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s) |
---|
| 41 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s) |
---|
| 42 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
---|
| 43 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
---|
| 44 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
---|
| 48 | |
---|
| 49 | REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: Km ! Exchange coefficient for momentum at interface between layers |
---|
| 50 | REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: Kh ! Exchange coefficient for heat flux at interface between layers |
---|
| 51 | |
---|
| 52 | ! Local variables |
---|
| 53 | REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange |
---|
| 54 | REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface |
---|
| 55 | REAL, DIMENSION(ngrid,nlay+1) :: dz_interf ! distance between two consecutive interfaces |
---|
| 56 | REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number |
---|
| 57 | REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number |
---|
| 58 | REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum |
---|
| 59 | REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat |
---|
| 60 | |
---|
| 61 | INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) |
---|
| 62 | REAL :: Ri0 ! parameter for Sm stability function |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | ! Initializations: |
---|
| 67 | !================ |
---|
| 68 | |
---|
| 69 | DO igrid=1,ngrid |
---|
| 70 | z_interf(igrid,1) = 0.0 |
---|
| 71 | Ri(igrid,1) = 0.1! TO DO |
---|
| 72 | Sm(igrid,1) = 0.0 ! TO DO |
---|
| 73 | END DO |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | ! Calculation of altitude of layers' bottom interfaces: |
---|
| 78 | !===================================================== |
---|
| 79 | |
---|
| 80 | DO ilay=2,nlay+1 |
---|
| 81 | DO igrid=1,ngrid |
---|
| 82 | dz_interf(igrid,ilay) = rg * temp(igrid,ilay) / rg * log(play(igrid,ilay-1)/play(igrid,ilay)) |
---|
| 83 | z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay) |
---|
| 84 | ENDDO |
---|
| 85 | ENDDO |
---|
| 86 | |
---|
| 87 | |
---|
| 88 | ! Computing the mixing length: |
---|
| 89 | !============================= |
---|
| 90 | |
---|
| 91 | DO ilay=2,nlay+1 |
---|
| 92 | DO igrid=1,ngrid |
---|
| 93 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
| 94 | ENDDO |
---|
| 95 | ENDDO |
---|
| 96 | |
---|
| 97 | ! Computes the gradient Richardson's number and stability functions: |
---|
| 98 | !=================================================================== |
---|
| 99 | |
---|
| 100 | DO ilay=2,nlay+1 |
---|
| 101 | DO igrid=1,ngrid |
---|
| 102 | ! Warning: sure that gradient of temp and wind should be calculated with dz_interf and not with dz_lay? |
---|
| 103 | Ri(igrid,ilay) = rg * (temp(igrid,ilay) - temp(igrid,ilay-1)) / dz_interf(igrid,ilay) / & |
---|
| 104 | (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 + & |
---|
| 105 | ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 ) |
---|
| 106 | |
---|
| 107 | ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 |
---|
| 108 | Ri0=2.*ric*cinf/rpi/cn |
---|
| 109 | |
---|
| 110 | IF (Ri(igrid,ilay) < 0.) THEN |
---|
| 111 | Sm(igrid,ilay) = 2./rpi * cinf * atan(-Ri(igrid,ilay)/Ri0) + cn |
---|
| 112 | Prandtl(igrid,ilay) = 1. + Ri(igrid,ilay) * pr_slope |
---|
| 113 | ELSE |
---|
| 114 | Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric)) |
---|
| 115 | Prandtl(igrid,ilay) = 2./rpi * atan(-Ri(igrid,ilay)) - 1. + pr_asym ! exp(Ri(igrid,ilay)) |
---|
| 116 | END IF |
---|
| 117 | |
---|
| 118 | Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay) |
---|
| 119 | |
---|
| 120 | ENDDO |
---|
| 121 | ENDDO |
---|
| 122 | |
---|
| 123 | ! Computing the TKE: |
---|
| 124 | !=================== |
---|
| 125 | IF (iflag_atke == 0) THEN |
---|
| 126 | |
---|
| 127 | ! stationary solution neglecting the vertical transport of TKE by turbulence |
---|
| 128 | DO ilay=2,nlay+1 |
---|
| 129 | DO igrid=1,ngrid |
---|
| 130 | tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & |
---|
| 131 | (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 + & |
---|
| 132 | ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 ) * & |
---|
| 133 | (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) |
---|
| 134 | ENDDO |
---|
| 135 | ENDDO |
---|
| 136 | |
---|
| 137 | ELSE ! TO DO |
---|
| 138 | |
---|
[4463] | 139 | call abort_physic("atke_compute_km_kh", & |
---|
| 140 | 'traitement non-stationnaire de la tke pas encore prevu', 1) |
---|
[4449] | 141 | |
---|
| 142 | END IF |
---|
| 143 | |
---|
| 144 | |
---|
| 145 | ! Computing eddy diffusivity coefficients: |
---|
| 146 | !======================================== |
---|
| 147 | DO ilay=2,nlay+1 |
---|
| 148 | DO igrid=1,ngrid |
---|
| 149 | Km(igrid,ilay) = l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
| 150 | Kh(igrid,ilay) = l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
| 151 | END DO |
---|
| 152 | END DO |
---|
| 153 | |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | end subroutine atke_compute_km_kh |
---|
| 157 | |
---|
| 158 | |
---|
| 159 | end module atke_exchange_coeff_mod |
---|