| 1 | !======================================================================= |
|---|
| 2 | ! CALLTHERM_INTERFACE |
|---|
| 3 | !======================================================================= |
|---|
| 4 | ! Main interface to the Martian thermal plume model |
|---|
| 5 | ! This interface handles sub-timesteps for this model |
|---|
| 6 | ! A call to this interface must be inserted in the main 'physics' routine |
|---|
| 7 | ! NB: for information: |
|---|
| 8 | ! In the Mars LMD-GCM, the thermal plume model is called after |
|---|
| 9 | ! the vertical turbulent mixing scheme (Mellor and Yamada) |
|---|
| 10 | ! and the surface layer scheme (Richardson-based surface layer + subgrid gustiness) |
|---|
| 11 | ! Other routines called before the thermals model are |
|---|
| 12 | ! radiative transfer and (orographic) gravity wave drag. |
|---|
| 13 | ! ----------------------------------------------------------------------- |
|---|
| 14 | ! Author : A. Colaitis 2011-01-05 (with updates 2011-2013) |
|---|
| 15 | ! after C. Rio and F. Hourdin |
|---|
| 16 | ! Institution : Laboratoire de Meteorologie Dynamique (LMD) Paris, France |
|---|
| 17 | ! ----------------------------------------------------------------------- |
|---|
| 18 | ! Corresponding author : A. Spiga aymeric.spiga_AT_upmc.fr |
|---|
| 19 | ! ----------------------------------------------------------------------- |
|---|
| 20 | ! ASSOCIATED FILES |
|---|
| 21 | ! --> thermcell_main_mars.F90 |
|---|
| 22 | ! --> thermcell_dqup.F90 |
|---|
| 23 | ! --> comtherm_h.F90 |
|---|
| 24 | ! ----------------------------------------------------------------------- |
|---|
| 25 | ! Reference paper: |
|---|
| 26 | ! A. Colaïtis, A. Spiga, F. Hourdin, C. Rio, F. Forget, and E. Millour. |
|---|
| 27 | ! A thermal plume model for the Martian convective boundary layer. |
|---|
| 28 | ! Journal of Geophysical Research (Planets), 118:1468-1487, July 2013. |
|---|
| 29 | ! http://dx.doi.org/10.1002/jgre.20104 |
|---|
| 30 | ! http://arxiv.org/abs/1306.6215 |
|---|
| 31 | ! ----------------------------------------------------------------------- |
|---|
| 32 | ! Reference paper for terrestrial plume model: |
|---|
| 33 | ! C. Rio and F. Hourdin. |
|---|
| 34 | ! A thermal plume model for the convective boundary layer : Representation of cumulus clouds. |
|---|
| 35 | ! Journal of the Atmospheric Sciences, 65:407-425, 2008. |
|---|
| 36 | ! ----------------------------------------------------------------------- |
|---|
| 37 | |
|---|
| 38 | SUBROUTINE calltherm_interface (ngrid,nlayer,nq, & |
|---|
| 39 | & tracer,igcm_co2, & |
|---|
| 40 | & zzlev,zzlay, & |
|---|
| 41 | & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, & |
|---|
| 42 | & pplay,pplev,pphi,zpopsk, & |
|---|
| 43 | & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke, & |
|---|
| 44 | & pdhdif,hfmax,wstar,sensibFlux) |
|---|
| 45 | |
|---|
| 46 | use comtherm_h |
|---|
| 47 | use tracer_mod, only: nqmx,noms |
|---|
| 48 | |
|---|
| 49 | ! SHARED VARIABLES. This needs adaptations in another climate model. |
|---|
| 50 | ! contains physical constant values such as |
|---|
| 51 | ! "g" : gravitational acceleration (m.s-2) |
|---|
| 52 | ! "r" : recuced gas constant (J.K-1.mol-1) |
|---|
| 53 | ! "cpp" : specific heat of the atmosphere (J.kg-1.K-1) |
|---|
| 54 | USE comcstfi_h |
|---|
| 55 | |
|---|
| 56 | implicit none |
|---|
| 57 | |
|---|
| 58 | !-------------------------------------------------------- |
|---|
| 59 | ! Input Variables |
|---|
| 60 | !-------------------------------------------------------- |
|---|
| 61 | |
|---|
| 62 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points |
|---|
| 63 | INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points |
|---|
| 64 | INTEGER, INTENT(IN) :: nq ! number of tracer species |
|---|
| 65 | REAL, INTENT(IN) :: ptimestep !timestep (s) |
|---|
| 66 | REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) !intermediate pressure levels (Pa) |
|---|
| 67 | REAL, INTENT(IN) :: pplay(ngrid,nlayer) !Pressure at the middle of the layers (Pa) |
|---|
| 68 | REAL, INTENT(IN) :: pphi(ngrid,nlayer) !Geopotential at the middle of the layers (m2s-2) |
|---|
| 69 | REAL, INTENT(IN) :: pu(ngrid,nlayer),pv(ngrid,nlayer) !u,v components of the wind (ms-1) |
|---|
| 70 | REAL, INTENT(IN) :: pt(ngrid,nlayer),pq(ngrid,nlayer,nq)!temperature (K) and tracer concentration (kg/kg) |
|---|
| 71 | REAL, INTENT(IN) :: zzlay(ngrid,nlayer) ! altitude at the middle of the layers |
|---|
| 72 | REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries |
|---|
| 73 | LOGICAL, INTENT(IN) :: tracer ! =.true. if tracers are present and to be transported |
|---|
| 74 | INTEGER, INTENT(IN) :: igcm_co2 ! index of the CO2 tracer in mixing ratio array |
|---|
| 75 | ! --> 0 if no tracer is CO2 (or no tracer at all) |
|---|
| 76 | ! --> this prepares special treatment for polar night mixing |
|---|
| 77 | ! (see thermcell_main_mars) |
|---|
| 78 | REAL, INTENT(IN) :: pdu(ngrid,nlayer),pdv(ngrid,nlayer) ! wind velocity change from routines called |
|---|
| 79 | ! before thermals du/dt (m/s/s) |
|---|
| 80 | REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tracer concentration change from routines called |
|---|
| 81 | ! before thermals dq/dt (kg/kg/s) |
|---|
| 82 | REAL, INTENT(IN) :: pdt(ngrid,nlayer) ! temperature change from routines called before thermals dT/dt (K/s) |
|---|
| 83 | REAL, INTENT(IN) :: q2(ngrid,nlayer+1) ! turbulent kinetic energy |
|---|
| 84 | REAL, INTENT(IN) :: zpopsk(ngrid,nlayer) ! ratio of pressure at middle of layer to surface pressure, |
|---|
| 85 | ! to the power r/cp, i.e. zpopsk=(pplay(ig,l)/pplev(ig,1))**rcp |
|---|
| 86 | REAL, INTENT(IN) :: pdhdif(ngrid,nlayer) ! potential temperature change from turbulent diffusion scheme dT/dt (K/s) |
|---|
| 87 | REAL, INTENT(IN) :: sensibFlux(ngrid) ! sensible heat flux computed from surface layer scheme |
|---|
| 88 | |
|---|
| 89 | !-------------------------------------------------------- |
|---|
| 90 | ! Output Variables |
|---|
| 91 | !-------------------------------------------------------- |
|---|
| 92 | |
|---|
| 93 | REAL, INTENT(OUT) :: pdu_th(ngrid,nlayer) ! wind velocity change from thermals du/dt (m/s/s) |
|---|
| 94 | REAL, INTENT(OUT) :: pdv_th(ngrid,nlayer) ! wind velocity change from thermals dv/dt (m/s/s) |
|---|
| 95 | REAL, INTENT(OUT) :: pdt_th(ngrid,nlayer) ! temperature change from thermals dT/dt (K/s) |
|---|
| 96 | REAL, INTENT(OUT) :: pdq_th(ngrid,nlayer,nq) ! tracer change from thermals dq/dt (kg/kg/s) |
|---|
| 97 | INTEGER, INTENT(OUT) :: lmax(ngrid) ! layer number reacher by thermals in grid point |
|---|
| 98 | REAL, INTENT(OUT) :: zmaxth(ngrid) ! equivalent to lmax, but in (m), interpolated |
|---|
| 99 | REAL, INTENT(OUT) :: pbl_dtke(ngrid,nlayer+1) ! turbulent kinetic energy change from thermals dtke/dt |
|---|
| 100 | REAL, INTENT(OUT) :: wstar(ngrid) ! free convection velocity (m/s) |
|---|
| 101 | |
|---|
| 102 | !-------------------------------------------------------- |
|---|
| 103 | ! Thermals local variables |
|---|
| 104 | !-------------------------------------------------------- |
|---|
| 105 | REAL zu(ngrid,nlayer), zv(ngrid,nlayer) |
|---|
| 106 | REAL zt(ngrid,nlayer) |
|---|
| 107 | REAL d_t_ajs(ngrid,nlayer) |
|---|
| 108 | REAL d_u_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq) |
|---|
| 109 | REAL d_v_ajs(ngrid,nlayer) |
|---|
| 110 | REAL fm_therm(ngrid,nlayer+1), entr_therm(ngrid,nlayer) |
|---|
| 111 | REAL detr_therm(ngrid,nlayer),detrmod(ngrid,nlayer) |
|---|
| 112 | REAL zw2(ngrid,nlayer+1) |
|---|
| 113 | REAL fraca(ngrid,nlayer+1),zfraca(ngrid,nlayer+1) |
|---|
| 114 | REAL q_therm(ngrid,nlayer), pq_therm(ngrid,nlayer,nq) |
|---|
| 115 | REAL q2_therm(ngrid,nlayer), dq2_therm(ngrid,nlayer) |
|---|
| 116 | REAL lmax_real(ngrid) |
|---|
| 117 | REAL masse(ngrid,nlayer) |
|---|
| 118 | |
|---|
| 119 | INTEGER l,ig,iq,ii(1),k |
|---|
| 120 | CHARACTER (LEN=20) modname |
|---|
| 121 | |
|---|
| 122 | !-------------------------------------------------------- |
|---|
| 123 | ! Local variables for sub-timestep |
|---|
| 124 | !-------------------------------------------------------- |
|---|
| 125 | |
|---|
| 126 | REAL d_t_the(ngrid,nlayer), d_q_the(ngrid,nlayer,nq) |
|---|
| 127 | INTEGER isplit |
|---|
| 128 | REAL fact |
|---|
| 129 | REAL zfm_therm(ngrid,nlayer+1),zdt |
|---|
| 130 | REAL zentr_therm(ngrid,nlayer),zdetr_therm(ngrid,nlayer) |
|---|
| 131 | REAL zheatFlux(ngrid,nlayer) |
|---|
| 132 | REAL zheatFlux_down(ngrid,nlayer) |
|---|
| 133 | REAL zbuoyancyOut(ngrid,nlayer) |
|---|
| 134 | REAL zbuoyancyEst(ngrid,nlayer) |
|---|
| 135 | REAL zzw2(ngrid,nlayer+1) |
|---|
| 136 | REAL zmax(ngrid) |
|---|
| 137 | INTEGER ndt,limz |
|---|
| 138 | |
|---|
| 139 | !-------------------------------------------------------- |
|---|
| 140 | ! Diagnostics |
|---|
| 141 | !-------------------------------------------------------- |
|---|
| 142 | |
|---|
| 143 | REAL heatFlux(ngrid,nlayer) |
|---|
| 144 | REAL heatFlux_down(ngrid,nlayer) |
|---|
| 145 | REAL buoyancyOut(ngrid,nlayer) |
|---|
| 146 | REAL buoyancyEst(ngrid,nlayer) |
|---|
| 147 | REAL hfmax(ngrid),wmax(ngrid) |
|---|
| 148 | REAL pbl_teta(ngrid),dteta(ngrid,nlayer) |
|---|
| 149 | REAL rpdhd(ngrid,nlayer) |
|---|
| 150 | REAL wtdif(ngrid,nlayer),rho(ngrid,nlayer) |
|---|
| 151 | REAL wtth(ngrid,nlayer) |
|---|
| 152 | |
|---|
| 153 | ! ********************************************************************** |
|---|
| 154 | ! Initializations |
|---|
| 155 | ! ********************************************************************** |
|---|
| 156 | |
|---|
| 157 | lmax(:)=0 |
|---|
| 158 | pdu_th(:,:)=0. |
|---|
| 159 | pdv_th(:,:)=0. |
|---|
| 160 | pdt_th(:,:)=0. |
|---|
| 161 | entr_therm(:,:)=0. |
|---|
| 162 | detr_therm(:,:)=0. |
|---|
| 163 | q2_therm(:,:)=0. |
|---|
| 164 | dq2_therm(:,:)=0. |
|---|
| 165 | pbl_dtke(:,:)=0. |
|---|
| 166 | fm_therm(:,:)=0. |
|---|
| 167 | zw2(:,:)=0. |
|---|
| 168 | fraca(:,:)=0. |
|---|
| 169 | zfraca(:,:)=0. |
|---|
| 170 | if (tracer) then |
|---|
| 171 | pdq_th(:,:,:)=0. |
|---|
| 172 | end if |
|---|
| 173 | d_t_ajs(:,:)=0. |
|---|
| 174 | d_u_ajs(:,:)=0. |
|---|
| 175 | d_v_ajs(:,:)=0. |
|---|
| 176 | d_q_ajs(:,:,:)=0. |
|---|
| 177 | heatFlux(:,:)=0. |
|---|
| 178 | heatFlux_down(:,:)=0. |
|---|
| 179 | buoyancyOut(:,:)=0. |
|---|
| 180 | buoyancyEst(:,:)=0. |
|---|
| 181 | zmaxth(:)=0. |
|---|
| 182 | lmax_real(:)=0. |
|---|
| 183 | |
|---|
| 184 | |
|---|
| 185 | ! ********************************************************************** |
|---|
| 186 | ! Preparing inputs for the thermals: increment tendancies |
|---|
| 187 | ! from other subroutines called before the thermals model |
|---|
| 188 | ! ********************************************************************** |
|---|
| 189 | |
|---|
| 190 | zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep ! u-component of wind velocity |
|---|
| 191 | zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep ! v-component of wind velocity |
|---|
| 192 | zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep ! temperature |
|---|
| 193 | |
|---|
| 194 | pq_therm(:,:,:)=0. ! tracer concentration |
|---|
| 195 | |
|---|
| 196 | if(qtransport_thermals) then |
|---|
| 197 | if(tracer) then ! tracer is a logical that is true if tracer must be transported in the GCM physics |
|---|
| 198 | pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep ! tracer concentration |
|---|
| 199 | endif |
|---|
| 200 | endif |
|---|
| 201 | |
|---|
| 202 | |
|---|
| 203 | IF(dtke_thermals) THEN |
|---|
| 204 | DO l=1,nlayer |
|---|
| 205 | q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1)) |
|---|
| 206 | ENDDO |
|---|
| 207 | ENDIF |
|---|
| 208 | |
|---|
| 209 | ! ********************************************************************** |
|---|
| 210 | ! --> CALLTHERM |
|---|
| 211 | ! SUB-TIMESTEP LOOP |
|---|
| 212 | ! ********************************************************************** |
|---|
| 213 | |
|---|
| 214 | zdt=ptimestep/REAL(nsplit_thermals) !subtimestep |
|---|
| 215 | |
|---|
| 216 | DO isplit=1,nsplit_thermals !temporal loop on the subtimestep |
|---|
| 217 | |
|---|
| 218 | ! Initialization of intermediary variables |
|---|
| 219 | |
|---|
| 220 | zzw2(:,:)=0. |
|---|
| 221 | zmax(:)=0. |
|---|
| 222 | lmax(:)=0 |
|---|
| 223 | |
|---|
| 224 | if (nq .ne. 0 .and. igcm_co2 .ne. 0) then !initialize co2 tracer tendancy |
|---|
| 225 | d_q_the(:,:,igcm_co2)=0. |
|---|
| 226 | endif |
|---|
| 227 | |
|---|
| 228 | ! CALL to main thermal routine |
|---|
| 229 | CALL thermcell_main_mars(ngrid,nlayer,nq & |
|---|
| 230 | & ,tracer,igcm_co2 & |
|---|
| 231 | & ,zdt & |
|---|
| 232 | & ,pplay,pplev,pphi,zzlev,zzlay & |
|---|
| 233 | & ,zu,zv,zt,pq_therm,q2_therm & |
|---|
| 234 | & ,d_t_the,d_q_the & |
|---|
| 235 | & ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax,limz & |
|---|
| 236 | & ,zzw2,fraca,zpopsk & |
|---|
| 237 | & ,zheatFlux,zheatFlux_down & |
|---|
| 238 | & ,zbuoyancyOut,zbuoyancyEst) |
|---|
| 239 | |
|---|
| 240 | fact=1./REAL(nsplit_thermals) |
|---|
| 241 | |
|---|
| 242 | ! Update thermals tendancies |
|---|
| 243 | d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact !temperature |
|---|
| 244 | if (igcm_co2 .ne. 0) then |
|---|
| 245 | d_q_the(:,:,igcm_co2)=d_q_the(:,:,igcm_co2)*ptimestep*fact !co2 mass mixing ratio |
|---|
| 246 | endif |
|---|
| 247 | zmaxth(:)=zmaxth(:)+zmax(:)*fact !thermals height |
|---|
| 248 | lmax_real(:)=lmax_real(:)+float(lmax(:))*fact !thermals height index |
|---|
| 249 | fm_therm(:,:)=fm_therm(:,:) & !upward mass flux |
|---|
| 250 | & +zfm_therm(:,:)*fact |
|---|
| 251 | entr_therm(:,:)=entr_therm(:,:) & !entrainment mass flux |
|---|
| 252 | & +zentr_therm(:,:)*fact |
|---|
| 253 | detr_therm(:,:)=detr_therm(:,:) & !detrainment mass flux |
|---|
| 254 | & +zdetr_therm(:,:)*fact |
|---|
| 255 | zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact !updraft fractional coverage |
|---|
| 256 | heatFlux(:,:)=heatFlux(:,:) & !upward heat flux |
|---|
| 257 | & +zheatFlux(:,:)*fact |
|---|
| 258 | heatFlux_down(:,:)=heatFlux_down(:,:) & !downward heat flux |
|---|
| 259 | & +zheatFlux_down(:,:)*fact |
|---|
| 260 | buoyancyOut(:,:)=buoyancyOut(:,:) & !plume final buoyancy |
|---|
| 261 | & +zbuoyancyOut(:,:)*fact |
|---|
| 262 | buoyancyEst(:,:)=buoyancyEst(:,:) & !plume estimated buoyancy used for vertical velocity computation |
|---|
| 263 | & +zbuoyancyEst(:,:)*fact |
|---|
| 264 | zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact !vertical velocity |
|---|
| 265 | |
|---|
| 266 | ! Save tendancies |
|---|
| 267 | d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) !temperature tendancy (delta T) |
|---|
| 268 | if (igcm_co2 .ne. 0) then |
|---|
| 269 | d_q_ajs(:,:,igcm_co2)=d_q_ajs(:,:,igcm_co2)+d_q_the(:,:,igcm_co2) !tracer tendancy (delta q) |
|---|
| 270 | endif |
|---|
| 271 | |
|---|
| 272 | ! Increment temperature and co2 concentration for next pass in subtimestep loop |
|---|
| 273 | zt(:,:) = zt(:,:) + d_t_the(:,:) !temperature |
|---|
| 274 | if (igcm_co2 .ne. 0) then |
|---|
| 275 | pq_therm(:,:,igcm_co2) = & |
|---|
| 276 | & pq_therm(:,:,igcm_co2) + d_q_the(:,:,igcm_co2) !co2 tracer |
|---|
| 277 | endif |
|---|
| 278 | |
|---|
| 279 | |
|---|
| 280 | ENDDO ! isplit |
|---|
| 281 | !**************************************************************** |
|---|
| 282 | |
|---|
| 283 | ! Now that we have computed total entrainment and detrainment, we can |
|---|
| 284 | ! advect u, v, and q in thermals. (potential temperature and co2 MMR |
|---|
| 285 | ! have already been advected in thermcell_main because they are coupled |
|---|
| 286 | ! to the determination of the thermals caracteristics). This is done |
|---|
| 287 | ! separatly because u,v, and q are not used in thermcell_main for |
|---|
| 288 | ! any thermals-related computation : they are purely passive. |
|---|
| 289 | |
|---|
| 290 | ! mass of cells |
|---|
| 291 | do l=1,nlayer |
|---|
| 292 | masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g |
|---|
| 293 | enddo |
|---|
| 294 | |
|---|
| 295 | ! recompute detrainment mass flux from entrainment and upward mass flux |
|---|
| 296 | ! this ensure mass flux conservation |
|---|
| 297 | detrmod(:,:)=0. |
|---|
| 298 | do l=1,limz |
|---|
| 299 | do ig=1,ngrid |
|---|
| 300 | detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) & |
|---|
| 301 | & +entr_therm(ig,l) |
|---|
| 302 | if (detrmod(ig,l).lt.0.) then |
|---|
| 303 | entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l) |
|---|
| 304 | detrmod(ig,l)=0. |
|---|
| 305 | endif |
|---|
| 306 | enddo |
|---|
| 307 | enddo |
|---|
| 308 | |
|---|
| 309 | ! u component of wind velocity advection in thermals |
|---|
| 310 | ! result is a derivative (d_u_ajs in m/s/s) |
|---|
| 311 | ndt=10 |
|---|
| 312 | call thermcell_dqup(ngrid,nlayer,ptimestep & |
|---|
| 313 | & ,fm_therm,entr_therm,detrmod, & |
|---|
| 314 | & masse,zu,d_u_ajs,ndt,limz) |
|---|
| 315 | |
|---|
| 316 | ! v component of wind velocity advection in thermals |
|---|
| 317 | ! result is a derivative (d_v_ajs in m/s/s) |
|---|
| 318 | call thermcell_dqup(ngrid,nlayer,ptimestep & |
|---|
| 319 | & ,fm_therm,entr_therm,detrmod, & |
|---|
| 320 | & masse,zv,d_v_ajs,ndt,limz) |
|---|
| 321 | |
|---|
| 322 | ! non co2 tracers advection in thermals |
|---|
| 323 | ! result is a derivative (d_q_ajs in kg/kg/s) |
|---|
| 324 | |
|---|
| 325 | if (nq .ne. 0.) then |
|---|
| 326 | DO iq=1,nq |
|---|
| 327 | if (iq .ne. igcm_co2) then |
|---|
| 328 | call thermcell_dqup(ngrid,nlayer,ptimestep & |
|---|
| 329 | & ,fm_therm,entr_therm,detrmod, & |
|---|
| 330 | & masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),ndt,limz) |
|---|
| 331 | endif |
|---|
| 332 | ENDDO |
|---|
| 333 | endif |
|---|
| 334 | |
|---|
| 335 | ! tke advection in thermals |
|---|
| 336 | ! result is a tendancy (d_u_ajs in J) |
|---|
| 337 | if (dtke_thermals) then |
|---|
| 338 | call thermcell_dqup(ngrid,nlayer,ptimestep & |
|---|
| 339 | & ,fm_therm,entr_therm,detrmod, & |
|---|
| 340 | & masse,q2_therm,dq2_therm,ndt,limz) |
|---|
| 341 | endif |
|---|
| 342 | |
|---|
| 343 | ! compute wmax for diagnostics |
|---|
| 344 | DO ig=1,ngrid |
|---|
| 345 | wmax(ig)=MAXVAL(zw2(ig,:)) |
|---|
| 346 | ENDDO |
|---|
| 347 | |
|---|
| 348 | ! ********************************************************************** |
|---|
| 349 | ! ********************************************************************** |
|---|
| 350 | ! ********************************************************************** |
|---|
| 351 | ! CALLTHERM END |
|---|
| 352 | ! ********************************************************************** |
|---|
| 353 | ! ********************************************************************** |
|---|
| 354 | ! ********************************************************************** |
|---|
| 355 | |
|---|
| 356 | |
|---|
| 357 | ! ********************************************************************** |
|---|
| 358 | ! Preparing outputs |
|---|
| 359 | ! ********************************************************************** |
|---|
| 360 | |
|---|
| 361 | do l=1,limz |
|---|
| 362 | pdu_th(:,l)=d_u_ajs(:,l) |
|---|
| 363 | pdv_th(:,l)=d_v_ajs(:,l) |
|---|
| 364 | enddo |
|---|
| 365 | |
|---|
| 366 | ! if tracers are transported in thermals, update output variables, else these are 0. |
|---|
| 367 | if(qtransport_thermals) then |
|---|
| 368 | if(tracer) then |
|---|
| 369 | do iq=1,nq |
|---|
| 370 | if (iq .ne. igcm_co2) then |
|---|
| 371 | do l=1,limz |
|---|
| 372 | pdq_th(:,l,iq)=d_q_ajs(:,l,iq) !non-co2 tracers d_q_ajs are dq/dt (kg/kg/s) |
|---|
| 373 | enddo |
|---|
| 374 | else |
|---|
| 375 | do l=1,limz |
|---|
| 376 | pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep !co2 tracer d_q_ajs is dq (kg/kg) |
|---|
| 377 | enddo |
|---|
| 378 | endif |
|---|
| 379 | enddo |
|---|
| 380 | endif |
|---|
| 381 | endif |
|---|
| 382 | |
|---|
| 383 | ! if tke is transported in thermals, update output variable, else this is 0. |
|---|
| 384 | IF(dtke_thermals) THEN |
|---|
| 385 | DO l=2,nlayer |
|---|
| 386 | pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l)) |
|---|
| 387 | ENDDO |
|---|
| 388 | |
|---|
| 389 | pbl_dtke(:,1)=0.5*dq2_therm(:,1) |
|---|
| 390 | pbl_dtke(:,nlayer+1)=0. |
|---|
| 391 | ENDIF |
|---|
| 392 | |
|---|
| 393 | ! update output variable for temperature. d_t_ajs is delta T in (K), pdt_th is dT/dt in (K/s) |
|---|
| 394 | do l=1,limz |
|---|
| 395 | pdt_th(:,l)=d_t_ajs(:,l)/ptimestep |
|---|
| 396 | enddo |
|---|
| 397 | |
|---|
| 398 | |
|---|
| 399 | ! ********************************************************************** |
|---|
| 400 | ! SURFACE LAYER INTERFACE |
|---|
| 401 | ! Compute the free convection velocity w* scale for surface layer gustiness |
|---|
| 402 | ! speed parameterization. The computed value of w* will be used at the next |
|---|
| 403 | ! timestep to modify surface-atmosphere exchange fluxes, because the surface |
|---|
| 404 | ! layer scheme and diffusion are called BEFORE the thermals. (outside of these |
|---|
| 405 | ! routines) |
|---|
| 406 | ! ********************************************************************** |
|---|
| 407 | |
|---|
| 408 | ! Potential temperature gradient |
|---|
| 409 | |
|---|
| 410 | dteta(:,nlayer)=0. |
|---|
| 411 | DO l=1,nlayer-1 |
|---|
| 412 | DO ig=1, ngrid |
|---|
| 413 | dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l)) & |
|---|
| 414 | & /(zzlay(ig,l+1)-zzlay(ig,l)) |
|---|
| 415 | ENDDO |
|---|
| 416 | ENDDO |
|---|
| 417 | |
|---|
| 418 | ! Computation of the PBL mixed layer temperature |
|---|
| 419 | |
|---|
| 420 | DO ig=1, ngrid |
|---|
| 421 | ii=MINLOC(abs(dteta(ig,1:lmax(ig)))) |
|---|
| 422 | pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1)) |
|---|
| 423 | ENDDO |
|---|
| 424 | |
|---|
| 425 | ! In order to have an accurate w*, we must add the heat flux from the |
|---|
| 426 | ! diffusion scheme to the computation of the maximum heat flux hfmax |
|---|
| 427 | ! Here pdhdif is the potential temperature change from the diffusion |
|---|
| 428 | ! scheme (Mellor and Yamada, see paper section 6, paragraph 57) |
|---|
| 429 | |
|---|
| 430 | ! compute rho as it is after the diffusion |
|---|
| 431 | |
|---|
| 432 | rho(:,:)=pplay(:,:) & |
|---|
| 433 | & /(r*(pt(:,:)+pdhdif(:,:)*zpopsk(:,:)*ptimestep)) |
|---|
| 434 | |
|---|
| 435 | ! integrate -rho*pdhdif |
|---|
| 436 | |
|---|
| 437 | rpdhd(:,:)=0. |
|---|
| 438 | |
|---|
| 439 | DO ig=1,ngrid |
|---|
| 440 | DO l=1,lmax(ig) |
|---|
| 441 | rpdhd(ig,l)=0. |
|---|
| 442 | DO k=1,l |
|---|
| 443 | rpdhd(ig,l)=rpdhd(ig,l)-rho(ig,k)*pdhdif(ig,k)* & |
|---|
| 444 | & (zzlev(ig,k+1)-zzlev(ig,k)) |
|---|
| 445 | ENDDO |
|---|
| 446 | rpdhd(ig,l)=rpdhd(ig,l)-sensibFlux(ig)/cpp |
|---|
| 447 | ENDDO |
|---|
| 448 | ENDDO |
|---|
| 449 | |
|---|
| 450 | ! compute w'theta' (vertical turbulent flux of temperature) from |
|---|
| 451 | ! the diffusion scheme |
|---|
| 452 | |
|---|
| 453 | wtdif(:,:)=rpdhd(:,:)/rho(:,:) |
|---|
| 454 | |
|---|
| 455 | ! Now we compute the contribution of the thermals to w'theta': |
|---|
| 456 | ! compute rho as it is after the thermals |
|---|
| 457 | |
|---|
| 458 | rho(:,:)=pplay(:,:)/(r*(zt(:,:))) |
|---|
| 459 | |
|---|
| 460 | ! integrate -rho*pdhdif |
|---|
| 461 | |
|---|
| 462 | DO ig=1,ngrid |
|---|
| 463 | DO l=1,lmax(ig) |
|---|
| 464 | rpdhd(ig,l)=0. |
|---|
| 465 | DO k=1,l |
|---|
| 466 | rpdhd(ig,l)=rpdhd(ig,l)-rho(ig,k)*(pdt_th(ig,k)/zpopsk(ig,k))* & |
|---|
| 467 | & (zzlev(ig,k+1)-zzlev(ig,k)) |
|---|
| 468 | ENDDO |
|---|
| 469 | rpdhd(ig,l)=rpdhd(ig,l)+ & |
|---|
| 470 | & rho(ig,1)*(heatFlux(ig,1)+heatFlux_down(ig,1)) |
|---|
| 471 | ENDDO |
|---|
| 472 | ENDDO |
|---|
| 473 | rpdhd(:,nlayer)=0. |
|---|
| 474 | |
|---|
| 475 | ! compute w'teta' from thermals |
|---|
| 476 | |
|---|
| 477 | wtth(:,:)=rpdhd(:,:)/rho(:,:) |
|---|
| 478 | |
|---|
| 479 | ! Add vertical turbulent heat fluxes from the thermals and the diffusion scheme |
|---|
| 480 | ! and compute the maximum |
|---|
| 481 | |
|---|
| 482 | DO ig=1,ngrid |
|---|
| 483 | hfmax(ig)=MAXVAL(wtth(ig,:)+wtdif(ig,:)) |
|---|
| 484 | ENDDO |
|---|
| 485 | |
|---|
| 486 | ! Finally we can compute the free convection velocity scale |
|---|
| 487 | ! We follow Spiga et. al 2010 (QJRMS) |
|---|
| 488 | ! ------------ |
|---|
| 489 | |
|---|
| 490 | DO ig=1, ngrid |
|---|
| 491 | IF (zmax(ig) .gt. 0.) THEN |
|---|
| 492 | wstar(ig)=(g*zmaxth(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.) |
|---|
| 493 | ELSE |
|---|
| 494 | wstar(ig)=0. |
|---|
| 495 | ENDIF |
|---|
| 496 | ENDDO |
|---|
| 497 | |
|---|
| 498 | END |
|---|
