| 1 |
#include "CPP_OPTIONS.h" |
| 2 |
|
| 3 |
SUBROUTINE MPDRIVER (CURTIME, DELT) |
| 4 |
IMPLICIT NONE |
| 5 |
C-- |
| 6 |
C-- SUBROUTINE PDRIVER (CURTIME) |
| 7 |
C-- |
| 8 |
C-- Purpose: stand-alone driver for MIT physical parametrization routines |
| 9 |
C-- Input : CURTIME : current time (in seconds) |
| 10 |
C-- grid-point model fields in common block: PHYGR1 |
| 11 |
C-- forcing fields in common blocks : LSMASK, FORFIX, FORCIN |
| 12 |
C-- Output : Diagnosed convective variables in common block: MITPHYGR2 |
| 13 |
C-- Other diagnosed variables in common block: MITPHYGR3 |
| 14 |
C-- Physical param. tendencies in common block: MITPHYTEN |
| 15 |
C-- Surface and upper boundary fluxes in common block: MITFLUXES |
| 16 |
C-- |
| 17 |
C |
| 18 |
C============================================================================= |
| 19 |
C |
| 20 |
C ******************************************************************* |
| 21 |
C * * |
| 22 |
C * MIT Physics package * |
| 23 |
C * * |
| 24 |
C * Version 1.0 * |
| 25 |
C * (July 2000) * |
| 26 |
C * * |
| 27 |
C ******************************************************************* |
| 28 |
C |
| 29 |
C |
| 30 |
C Created: 07/2000 Sandrine Bony & Kerry Emanuel |
| 31 |
C Modification by Olivier Pauluis (07-2001 --- now) |
| 32 |
C |
| 33 |
C -- Important: |
| 34 |
C ------------- |
| 35 |
C |
| 36 |
C * this version (1.0) of the physics package is intended to be used over |
| 37 |
C a tropical aqua-planet (no sea ice, no distinction between land and |
| 38 |
C ocean surfaces, no specific cloudiness associated with mid-latitude |
| 39 |
C weather systems, etc.) |
| 40 |
C ... but nothing precludes to complete it later on ! |
| 41 |
C |
| 42 |
C * for a non-equatorial version of the model, don't forget to specify |
| 43 |
C the latitude (RLAT) seen by radiation. |
| 44 |
C |
| 45 |
C * if you change the vertical resolution, change it accordingly |
| 46 |
C in the file ../inc/dimensions.h (this is for the radiation code) |
| 47 |
C and also in convect43b.F (NLEV=...) |
| 48 |
C (not very satisfying...will have to be improved) |
| 49 |
C |
| 50 |
C |
| 51 |
C -- Among important differences with Franco Molteni's physics: |
| 52 |
C ------------------------------------------------------------- |
| 53 |
C |
| 54 |
C - opposite convention used for numbering the vertical levels |
| 55 |
C (here, the 1st atmospheric level is the nearest to the surface) |
| 56 |
C |
| 57 |
C - QG1 (specific humidity) is in kg/kg |
| 58 |
C |
| 59 |
C - all the physics package is called atmospheric column by |
| 60 |
C atmospheric column (note that the radiation code would |
| 61 |
C easily accept to work with 3D fields) |
| 62 |
C |
| 63 |
C |
| 64 |
C -- If you are looking for some documentation about the parameterizations |
| 65 |
C used in this package: |
| 66 |
C -------------------------------------------------------------------------- |
| 67 |
C |
| 68 |
C Bony S and K A Emanuel, 2001: A parameterization of the cloudiness |
| 69 |
C associated with cumulus convection; Evaluation using TOGA COARE |
| 70 |
C data. J. Atmos. Sci., accepted. |
| 71 |
C |
| 72 |
C Emanuel K A, 1991: A scheme for representing cumulus convection in |
| 73 |
C large-scale models. J. Atmos. Sci., 48, 2313-2335. |
| 74 |
C |
| 75 |
C Emanuel K A and M Zivkovic-Rothman, 1999: Development and evaluation |
| 76 |
C of a convection scheme for use in climate models. J. Atmos. Sci., |
| 77 |
C 56, 1766-1782 (see also http://www-paoc.mit.edu/~emanuel/home.html) |
| 78 |
C |
| 79 |
C Fouquart Y and B Bonnel, 1980: Computation of solar heating of the |
| 80 |
C Earth's atmosphere: a new parameterization. Beitr. Phys. Atmos., |
| 81 |
C 53, 35-62. |
| 82 |
C |
| 83 |
C Morcrette J-J, 1991: Radiation and cloud radiative properties in the |
| 84 |
C European Centre for Medium-Range Weather Forecasts forecasting |
| 85 |
C system. J. Geophys. Res., 96, 9121-9132. |
| 86 |
C |
| 87 |
C-- Please, kindly report any problem or suggestion to the authors: |
| 88 |
C |
| 89 |
C bony@wind.mit.edu, emanuel@texmex.mit.edu |
| 90 |
C |
| 91 |
C============================================================================= |
| 92 |
|
| 93 |
C Resolution parameters |
| 94 |
C |
| 95 |
#include "atparam.h" |
| 96 |
#include "atparam1.h" |
| 97 |
#include "MITPHYS_PARAMS.h" |
| 98 |
C |
| 99 |
INTEGER NLON, NLAT, NLEV, NGP |
| 100 |
PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT ) |
| 101 |
C |
| 102 |
C Constants + functions of sigma and latitude |
| 103 |
C |
| 104 |
#include "Lev_def.h" |
| 105 |
C |
| 106 |
C Model variables, tendencies and fluxes on gaussian grid |
| 107 |
C |
| 108 |
#include "com_mitphysvar.h" |
| 109 |
C |
| 110 |
C Diagnostics from the convection scheme: |
| 111 |
C |
| 112 |
#include "com_convect.h" |
| 113 |
C |
| 114 |
C Assign thermodynamical constants: |
| 115 |
C |
| 116 |
#include "thermo.h" |
| 117 |
C |
| 118 |
C Surface forcing fields (time-inv. or functions of seasonal cycle) |
| 119 |
C |
| 120 |
#include "com_forcing1.h" |
| 121 |
|
| 122 |
c Misc: |
| 123 |
|
| 124 |
INTEGER I,J,K |
| 125 |
_RL CURTIME |
| 126 |
_RL DELT, ES, AUX, RSS, RS(NLEV), QMIN |
| 127 |
PARAMETER (QMIN=1.e-12) ! minimum humidity when negative |
| 128 |
C PARAMETER (DELT=600.D0) ! timestep in seconds (temporary: deltaT better) |
| 129 |
C (not very satisfying...will have to be improved) |
| 130 |
_RL TS(NGP) |
| 131 |
|
| 132 |
_RL SX(NLEV),HX(NLEV),LV(NLEV),GZ(NLEV),TVY,TVX,CPN(NLEV) |
| 133 |
|
| 134 |
c Convection: |
| 135 |
|
| 136 |
INTEGER IFLAG, NTRA, NL |
| 137 |
PARAMETER (NTRA=1) ! number of tracors in convection |
| 138 |
_RL CBMF, PRECIP, WD, TPRIME, QPRIME |
| 139 |
_RL TCONV(NLEV), RCONV(NLEV), RVCONV(NLEV) |
| 140 |
_RL TMEMO(NLEV), RMEMO(NLEV) |
| 141 |
_RL UCONV_E(NLEV), UCONV_W(NLEV) |
| 142 |
_RL VCONV_S(NLEV), VCONV_N(NLEV) |
| 143 |
_RL UCONV(NLEV), VCONV(NLEV) |
| 144 |
_RL TRACONV(NLEV,NTRA) |
| 145 |
_RL QCONDC(NLEV), FTCONV(NLEV), FRCONV(NLEV) |
| 146 |
_RL FUCONV(NLEV), FVCONV(NLEV), FTRACONV(NLEV,NTRA) |
| 147 |
C modif omp |
| 148 |
_RL PLCL |
| 149 |
|
| 150 |
c Soil moisture: bucket model |
| 151 |
|
| 152 |
_RL LANDMUL(NGP), FACTOR, EBKT, PBKT |
| 153 |
|
| 154 |
c Cloudiness and large-scale condensation: |
| 155 |
|
| 156 |
_RL PRADJ, CLDT, CLDWP |
| 157 |
_RL CLDF(NLEV), CLDQ(NLEV), FTADJ(NLEV), FRADJ(NLEV) |
| 158 |
_RL CLDEMI(NLEV), CLDTAU(NLEV), CLDFICE(NLEV) |
| 159 |
_RL CLDFRAD(NLEV), CLDEMIRAD(NLEV), CLDTAURAD(NLEV) |
| 160 |
|
| 161 |
|
| 162 |
_RL CLDT_M, CLDWP_M |
| 163 |
_RL CLDF_M(NLEV), CLDQ_M(NLEV) |
| 164 |
_RL CLDEMI_M(NLEV), CLDTAU_M(NLEV), CLDFICE_M(NLEV) |
| 165 |
|
| 166 |
|
| 167 |
c Radiation: |
| 168 |
|
| 169 |
c modif omp: now included in MITPHYS_PARAMS.h |
| 170 |
C _RL RADFREQ ! frequency of radiation calls |
| 171 |
C PARAMETER (RADFREQ=6*3600.) ! here: every 6 hours |
| 172 |
|
| 173 |
|
| 174 |
_RL SCON, CO2, ALB! solar cst, CO2 conc, sfc albedo |
| 175 |
PARAMETER (SCON=1370.0, CO2=330.0) !units: W/m2, ppm |
| 176 |
c PARAMETER (ALB=0.07) !units: fraction |
| 177 |
PARAMETER (ALB=0.10) !units: fraction |
| 178 |
|
| 179 |
C REPLACED by a logical DARAD in MITPHYS_PARAMS.h |
| 180 |
C CHARACTER*1 DARAD ! diurnally averaged radiation? |
| 181 |
C PARAMETER (DARAD='y') |
| 182 |
|
| 183 |
_RL RC0(NGP,NLEV) ! specified radiative heating rate |
| 184 |
_RL RCRF0(NLEV) ! specified net cld radiative forcing |
| 185 |
|
| 186 |
C modif omp: these parameters are defined in the MITPHYS namelist |
| 187 |
|
| 188 |
C LOGICAL RADINT, CRFINT, CLEARSKY |
| 189 |
C PARAMETER (RADINT=.TRUE.) ! TRUE=interactive radiation |
| 190 |
C PARAMETER (CRFINT=.TRUE.)! TRUE=interactive cld rad forcing |
| 191 |
C PARAMETER (CLEARSKY=.FALSE.)! TRUE=ignore cld-rad interactions |
| 192 |
C PARAMETER (RADINT=.FALSE.) ! TRUE=interactive radiation |
| 193 |
C PARAMETER (CRFINT=.FALSE.)! TRUE=interactive cld rad forcing |
| 194 |
C PARAMETER (CLEARSKY=.TRUE.)! TRUE=ignore cld-rad interactions |
| 195 |
|
| 196 |
INTEGER julien |
| 197 |
_RL rjulien |
| 198 |
|
| 199 |
_RL zlongi, dist, rmu0, fract, TSRAD, TIMEU, gmtime |
| 200 |
_RL radsol, albsol, albpla, RLON, RLAT, normalis |
| 201 |
_RL O3(NLEV), O3RAD(NLEV) |
| 202 |
_RL topsw, toplw, solsw, sollw |
| 203 |
_RL topsw0, toplw0, solsw0, sollw0, inso |
| 204 |
_RL heat(NLEV), heat0(NLEV), cool(NLEV), cool0(NLEV) |
| 205 |
_RL RADTIME |
| 206 |
SAVE RADTIME |
| 207 |
|
| 208 |
|
| 209 |
|
| 210 |
LOGICAL FirstCall |
| 211 |
DATA FirstCall /.TRUE./ |
| 212 |
SAVE FirstCall |
| 213 |
|
| 214 |
c Surface fluxes: |
| 215 |
|
| 216 |
C modif omp: these parameters are defined in the MITPHYS namelist |
| 217 |
|
| 218 |
C _RL CD ! drag coeff and sfc evap potential |
| 219 |
C PARAMETER (CD=1.2D-3) |
| 220 |
|
| 221 |
C _RL US0,VS0,WS0 ! mean and min sfc winds |
| 222 |
C PARAMETER (US0= 0.0,VS0=0.0,WS0=4.0) |
| 223 |
_RL EFRAC |
| 224 |
PARAMETER (EFRAC = 1.0) |
| 225 |
|
| 226 |
_RL US2,VS2,VSURF,VSPRIME,TC,ALV |
| 227 |
_RL TSA,PG,ROWS,DPH1I,FTSURF,FRSURF,SH,LH |
| 228 |
|
| 229 |
_RL SST_TEND, STL_TEND |
| 230 |
|
| 231 |
_RL DUMMY |
| 232 |
|
| 233 |
C modif omp: these parameters are defined in the MITPHYS namelist |
| 234 |
C OMP 07-18-01 |
| 235 |
C LOG_CORRECT: correction factor for wind distribution in the PBL. |
| 236 |
C Defined as the ratio of the 10m wind to the wind at the |
| 237 |
C lowest model level |
| 238 |
C LOG_CORRECT = 1: wind uniform in the lowest level |
| 239 |
C LOG_CORRECT = 2/3: logarithmic profile between the surface to 100m, with a |
| 240 |
C surfcae roughness ~ 10cm. |
| 241 |
|
| 242 |
C _RL LOG_CORRECT |
| 243 |
C PARAMETER ( LOG_CORRECT = 1.0 ) |
| 244 |
C PARAMETER ( LOG_CORRECT = 0.6666666666 ) |
| 245 |
|
| 246 |
|
| 247 |
|
| 248 |
C OMP: LOW LEVEL MIXING: |
| 249 |
|
| 250 |
C modif omp: these parameters are defined in the MITPHYS namelist |
| 251 |
C INTEGER K_PBL |
| 252 |
C PARAMETER (K_PBL = 8 ) |
| 253 |
C LOGICAL PBL_DIFF |
| 254 |
C PARAMETER (PBL_DIFF = .FALSE.) |
| 255 |
|
| 256 |
_RL U_PBLW, U_PBLE, V_PBLS, V_PBLN, PBL_INV |
| 257 |
_RL UW_FLUX(nlev), UE_FLUX(nlev), VN_FLUX(nlev),VS_FLUX(nlev) |
| 258 |
_RL DPI, VISC_PBL, DELTI_ADJ, DELTI_MIX |
| 259 |
_RL WIND_ADJ |
| 260 |
|
| 261 |
_RL FUSURF_W, FUSURF_E, FVSURF_S,FVSURF_N |
| 262 |
|
| 263 |
c Momentum drag: |
| 264 |
|
| 265 |
c _RL TAUM ! timescale for zonal wind relaxation |
| 266 |
c PARAMETER (TAUM=3.*86400.) ! 3 days |
| 267 |
|
| 268 |
_RL UMEAN, VMEAN ! specified mean barotropic flow |
| 269 |
PARAMETER (UMEAN=0.0, VMEAN=0.) ! for drag only |
| 270 |
|
| 271 |
C modif omp: these parameters are defined in the MITPHYS namelist |
| 272 |
C _RL CDM ! drag coeff for sfc flux of momentum |
| 273 |
C PARAMETER (CDM=1.2d-3) |
| 274 |
|
| 275 |
c _RL UZON(NLEV) |
| 276 |
|
| 277 |
C Atmospheric pressure levels used by the physics package: |
| 278 |
C PPLAY (mb): pressure surfaces where T, Q are computed |
| 279 |
C PAPRS (mb): pressure surfaces at the interface between layers |
| 280 |
C PPLAY1 and PAPRS1: same as PPLAY and PAPRS but in Pa |
| 281 |
C (should be improved later...use pSurfs instead) |
| 282 |
|
| 283 |
_RL DP, PPLAY(NLEV), PAPRS(NLEV+1) |
| 284 |
_RL PPLAY1(NLEV), PAPRS1(NLEV+1) |
| 285 |
|
| 286 |
c 5levels DATA PPLAY / 950.0, 775.0, 500.0, 250.0, 75.0 / |
| 287 |
c 5levels DATA PAPRS / 1000.0, 900.0, 650.0, 350.0, 150.0, 0.0 / |
| 288 |
c 5levels DATA O3 / 0.0, 0.0, 0.0, 0.0, 0.0 / |
| 289 |
|
| 290 |
DATA O3 / 0.048, 0.049, 0.051, 0.052, 0.053, 0.054, 0.056, |
| 291 |
: 0.056, 0.057, 0.058, 0.059, 0.059, 0.060, 0.060, |
| 292 |
: 0.061, 0.062, 0.063, 0.064, 0.066, 0.068, 0.069, |
| 293 |
: 0.072, 0.073, 0.075, 0.077, 0.081, 0.085, 0.090, |
| 294 |
: 0.096, 0.106, 0.117, 0.135, 0.153, 0.175, 0.202, |
| 295 |
: 0.243, 0.405, 1.000, 3.700, 10.124 / |
| 296 |
|
| 297 |
c mean CRF profile derived for SST=27C, U*=2m/s, U0=-5m/s, 360x40, 180E: |
| 298 |
DATA RCRF0 / 0.152041972, 0.183760509, 0.167694584, 0.16433309, |
| 299 |
: 0.16215764, 0.161808997, 0.163289174, 0.166284353, |
| 300 |
: 0.170786828, 0.176764682, 0.1845745, 0.19467622, |
| 301 |
: 0.201630697, 0.209879801, 0.223531589, 0.245427981, |
| 302 |
: 0.270941138, 0.301004857, 0.324049205, 0.340905219, |
| 303 |
: 0.360249668, 0.389321774, 0.412838638, 0.449643821, |
| 304 |
: 0.483938843, 0.525295079, 0.666303575, 0.935836792, |
| 305 |
: 1.25989616, 1.19505048, 1.40494013, 1.28129566, |
| 306 |
: -0.570507288, -0.167021349, -0.00548044546, 0.129761115, |
| 307 |
: -0.0882121176, -0.0795640424, -0.134504586, 0.0340689607 / |
| 308 |
|
| 309 |
C ========================================================================= |
| 310 |
C |
| 311 |
C Loop over atmospheric columns starts here: |
| 312 |
C |
| 313 |
C ========================================================================= |
| 314 |
|
| 315 |
CC (acz) check for SST |
| 316 |
CC write(6,*) 'beginning of mit_phydriver: SST1=',SST1(1) |
| 317 |
|
| 318 |
IF (FirstCall) THEN |
| 319 |
DO J = 1,NGP |
| 320 |
DO K = 1 , NLEV |
| 321 |
CLDF_MEAN(J,K) = 0.0 |
| 322 |
CLDQ_MEAN(J,K) = 0.0 |
| 323 |
END DO |
| 324 |
END DO |
| 325 |
END IF |
| 326 |
|
| 327 |
C note: using 1.0/DELT potentially unstable because of the |
| 328 |
C Adams-Bashford time-stepping. CAN be put to 1 if the physics is |
| 329 |
C taken out of the Adams-Bashford... |
| 330 |
|
| 331 |
|
| 332 |
IF (CONV_ADJ) THEN |
| 333 |
DELTI_ADJ = 1.0/ DELT |
| 334 |
WIND_ADJ = 0.6 |
| 335 |
ELSE |
| 336 |
DELTI_ADJ = 1.0/ DELT |
| 337 |
WIND_ADJ = 0.0 |
| 338 |
END IF |
| 339 |
DELTI_MIX = 0.6 / DELT |
| 340 |
|
| 341 |
|
| 342 |
|
| 343 |
RADTIME = RADTIME + DELT |
| 344 |
|
| 345 |
|
| 346 |
c -- pressure levels: |
| 347 |
|
| 348 |
DP = 1000./NLEV |
| 349 |
DO K = 1, NLEV |
| 350 |
PPLAY(K) = 1000. - (K-1)*DP |
| 351 |
IF(K.GT.1) PAPRS(K)=0.5*(PPLAY(K)+PPLAY(K-1)) |
| 352 |
ENDDO |
| 353 |
PAPRS(1) = 2.*PPLAY(1)-PAPRS(2) |
| 354 |
PAPRS(NLEV+1) = 2.*PAPRS(NLEV)-PAPRS(NLEV-1) |
| 355 |
PAPRS(NLEV+1) = MAX(0.1,PAPRS(NLEV+1)) |
| 356 |
|
| 357 |
c -- zonally-averaged zonal wind at each level: |
| 358 |
|
| 359 |
c DO K = 1, NLEV |
| 360 |
c UZON(K) = 0. |
| 361 |
c DO J = 1, NGP |
| 362 |
c UZON(K) = UZON(K) + UG1(J,K)/FLOAT(NGP) |
| 363 |
c ENDDO |
| 364 |
c ENDDO |
| 365 |
|
| 366 |
c -- loop over horizontal gridpoints: |
| 367 |
|
| 368 |
C Newtonian cooling rate |
| 369 |
|
| 370 |
IF (.NOT. RADINT ) THEN |
| 371 |
CALL NEWT_COOL(TG1,RC0) |
| 372 |
END IF |
| 373 |
|
| 374 |
|
| 375 |
DO 9999 J = 1, NGP |
| 376 |
|
| 377 |
c ----------------------------------------------------------------------- |
| 378 |
c 1. Preliminaries: |
| 379 |
c ----------------------------------------------------------------------- |
| 380 |
|
| 381 |
c -- check that humidities are not negative: |
| 382 |
DO K = 1, NLEV |
| 383 |
if (QG1(J,K).lt.0.0 ) then |
| 384 |
CC(acz) write(*,*) ' J, K, neg humidity: ',J, K,QG1(J,K) |
| 385 |
QG1(J,K) = MAX( QG1(J,K), QMIN ) |
| 386 |
endif |
| 387 |
ENDDO |
| 388 |
|
| 389 |
c -- initialize tendencies (MITPHYTEN common): |
| 390 |
|
| 391 |
DO K = 1, NLEV |
| 392 |
TT_CNV(J,K) = 0.0 |
| 393 |
QT_CNV(J,K) = 0.0 |
| 394 |
UT_CNV(J,K) = 0.0 |
| 395 |
VT_CNV(J,K) = 0.0 |
| 396 |
TT_LSC(J,K) = 0.0 |
| 397 |
QT_LSC(J,K) = 0.0 |
| 398 |
TT_PBL(J,K) = 0.0 |
| 399 |
QT_PBL(J,K) = 0.0 |
| 400 |
UT_PBL_E(J,K) = 0.0 |
| 401 |
UT_PBL_W(J,K) = 0.0 |
| 402 |
VT_PBL_S(J,K) = 0.0 |
| 403 |
VT_PBL_N(J,K) = 0.0 |
| 404 |
ENDDO |
| 405 |
PRECNV(J) = 0. |
| 406 |
PRECLS(J) = 0. |
| 407 |
LHF(J) = 0. |
| 408 |
SHF(J) = 0. |
| 409 |
|
| 410 |
c -- surface pressure (mb): |
| 411 |
|
| 412 |
PG = 0.01 * PSG1(J) |
| 413 |
|
| 414 |
c -- surface temperature (K): |
| 415 |
|
| 416 |
TS(J) =SST1(J)+FMASK1(J)*(STL1(J)-SST1(J)) |
| 417 |
|
| 418 |
CC (acz) check for SST |
| 419 |
CC write(6,*) 'start of phy_driver: J, FMASK1(J)=',J, FMASK1(J) |
| 420 |
CC if( J .EQ. NGP) THEN |
| 421 |
CC write(6,*) 'start of phy_driver: TS, SST1, STL1=',TS(NGP), SST1(NGP),STL1(NGP) |
| 422 |
CC end if |
| 423 |
|
| 424 |
c -- compute the land multiplicator factor for evaporation |
| 425 |
c (acz) adapted from (nprive) - Delworth and Manabe (1988) |
| 426 |
|
| 427 |
if (BUCKET(J) .GE. 0.75*BKT0) then |
| 428 |
FACTOR = 1 |
| 429 |
else |
| 430 |
FACTOR = BUCKET(J)/(0.75*BKT0) |
| 431 |
endif |
| 432 |
|
| 433 |
LANDMUL(J) = 1 + FMASK1(J) * (FACTOR-1) |
| 434 |
|
| 435 |
c -- compute saturation specific humidity: |
| 436 |
|
| 437 |
TC=TS(J)-273.15 |
| 438 |
ES=6.112*EXP(17.67*TC/(243.5+TC)) |
| 439 |
RSS=0.98*EPS*ES/(PG-(1.-EPS)*ES) |
| 440 |
|
| 441 |
DO K = 1, NLEV |
| 442 |
TC=TG1(J,K)-273.15 |
| 443 |
IF(TC.GE.0.0)THEN |
| 444 |
ES=6.112*EXP(17.67*TC/(243.5+TC)) |
| 445 |
ELSE |
| 446 |
ES=EXP(23.33086-6111.72784/TG1(J,K)+0.15215*LOG(TG1(J,K))) |
| 447 |
ENDIF |
| 448 |
ES=MIN(ES,PPLAY(K)) |
| 449 |
RS(K)=EPS*ES/(PPLAY(K)-(1.-EPS)*ES) |
| 450 |
ENDDO |
| 451 |
|
| 452 |
c -- eventually, prescribe the radiative cooling rate to a cst value: |
| 453 |
|
| 454 |
C omp: Prescribed cooling has been replaced by a newtonian cooling |
| 455 |
|
| 456 |
IF ( .NOT. RADINT ) THEN |
| 457 |
DO K = 1, NLEV |
| 458 |
C RC1(J,K) = - RC0(J,K)! net rad heating (units = K/day) |
| 459 |
RC1(J,K) = RC0(J,K) * 24.d0 * 3600.d0 ! net rad heating (units = K/day) |
| 460 |
ENDDO |
| 461 |
|
| 462 |
|
| 463 |
ENDIF |
| 464 |
|
| 465 |
c ----------------------------------------------------------------------- |
| 466 |
c 2. Convection subroutine |
| 467 |
c |
| 468 |
c CBMF : cloud-base mass flux (its value must be "remembered" by the |
| 469 |
c calling program between calls to the convection subroutine). |
| 470 |
c |
| 471 |
c QCONDC: in-cloud mixing ratio of condensed water produced by cumulus |
| 472 |
c convection [kg/kg], output of the convection scheme used in |
| 473 |
c the cloud parameterization. |
| 474 |
c |
| 475 |
c WD, TPRIME, QPRIME: convective downdraft velocity scale [m/s], |
| 476 |
c temperature perturbation scale [K] and specific humidity |
| 477 |
c perturbation scale [kg/kg] that will be used in surface |
| 478 |
c flux parameterizations. |
| 479 |
c ----------------------------------------------------------------------- |
| 480 |
|
| 481 |
c cloud-base mass flux: |
| 482 |
|
| 483 |
CBMF = CBMFG1(J) |
| 484 |
PLCL = 0 |
| 485 |
|
| 486 |
c NL: the maximum number of levels to which convection can penetrate, plus 1 |
| 487 |
c (NL MUST be less than or equal to NLEV-1) |
| 488 |
|
| 489 |
NL = 1 |
| 490 |
DO K = 1, NLEV |
| 491 |
IF(PPLAY(K).GE.78.0) NL = MAX(NL,K) |
| 492 |
ENDDO |
| 493 |
|
| 494 |
c T and Q profiles seen by the convection (and eventually altered by |
| 495 |
c the subroutine is dry convective adjustment occurs): |
| 496 |
|
| 497 |
DO K = 1, NLEV |
| 498 |
TCONV(K) = TG1(J,K) ! TCONV in K |
| 499 |
RCONV(K) = QG1(J,K) ! RCONV in kg/kg |
| 500 |
C OMP: Current implementation of the momentum transport consider the mean |
| 501 |
C velocity at the center of the cell. Future implementatioin may |
| 502 |
C separate between UGW, VGW, VGS, VGN. this would require modifying |
| 503 |
C convect43 to handle more than 2 velocity. |
| 504 |
|
| 505 |
UCONV(K) = 0.5 * (UGW(J,K) + UGE(J,K)) |
| 506 |
VCONV(K) = 0.5 * (VGS(J,K) + VGN(J,K)) |
| 507 |
|
| 508 |
UCONV_W(K) = UGW(J,K) |
| 509 |
UCONV_E(K) = UGE(J,K) |
| 510 |
VCONV_S(K) = VGS(J,K) |
| 511 |
VCONV_N(K) = VGN(J,K) |
| 512 |
DO I = 1, NTRA |
| 513 |
TRACONV(K,I) = 0.0 ! no tracors here |
| 514 |
ENDDO |
| 515 |
ENDDO |
| 516 |
|
| 517 |
c "official" version 4.3b of CONVECT plus diagnostics required for cld parameterization: |
| 518 |
|
| 519 |
|
| 520 |
C write(*,*) '--- convect: I2 = ', J, 'CBMF =', CBMF |
| 521 |
|
| 522 |
CALL CONVECT43C(TCONV, RCONV, RS, UCONV, VCONV, |
| 523 |
: TRACONV, |
| 524 |
: PPLAY, PAPRS, NLEV, NL, NTRA, DELT, |
| 525 |
: IFLAG, FTCONV, FRCONV, FUCONV, FVCONV, |
| 526 |
: FTRACONV, PRECIP, WD, TPRIME, QPRIME, QCONDC, |
| 527 |
: CBMF, PLCL ) |
| 528 |
C write(*,*) 'after convect, CBMF =', CBMF |
| 529 |
|
| 530 |
C if necessary, write to error file |
| 531 |
|
| 532 |
|
| 533 |
IF (IFLAG.EQ.4) THEN |
| 534 |
WRITE(*,120) CURTIME |
| 535 |
120 FORMAT(5X,'CFL condition violated at TIME ',e15.7) |
| 536 |
WRITE(*,*) 'LAT: ', LAT_G(J), 'LON: ', LON_G(J) |
| 537 |
ENDIF |
| 538 |
|
| 539 |
c T and Q tendencies associated with convection: |
| 540 |
|
| 541 |
|
| 542 |
|
| 543 |
DO K = 1, NL |
| 544 |
TT_CNV(J,K) = FTCONV(K) ! K/s |
| 545 |
QT_CNV(J,K) = FRCONV(K) ! kg/kg/s |
| 546 |
UT_CNV(J,K) = FUCONV(K) ! m/s^2 |
| 547 |
VT_CNV(J,K) = FVCONV(K) ! m/s^2 |
| 548 |
c |
| 549 |
c note that we do not update the wind, here |
| 550 |
c (should be done in principle if the dry adjustment |
| 551 |
c is called and if the wind-tendencies due to the |
| 552 |
c convection are to be taken into account). |
| 553 |
|
| 554 |
c modif omp: pass the convective adjustment as a time tendecy |
| 555 |
TT_ADJ(J,K) = DELTI_ADJ * (TCONV(K) - TG1(J,K)) |
| 556 |
QT_ADJ(J,K) = DELTI_ADJ * (RCONV(K) - QG1(J,K)) |
| 557 |
UT_ADJ(J,K) = WIND_ADJ * DELTI_ADJ * (UCONV(K) |
| 558 |
& - 0.5 * (UCONV_W(K) + UCONV_E(K))) |
| 559 |
VT_ADJ(J,K) = WIND_ADJ * DELTI_ADJ * (VCONV(K) |
| 560 |
& - 0.5 * (VCONV_S(K) + VCONV_N(K))) |
| 561 |
TG1(J,K) = TCONV(K) |
| 562 |
QG1(J,K) = RCONV(K) ! QG1 in kg/kg |
| 563 |
C UGW(J,K) = UCONV_W(K) |
| 564 |
C UGE(J,K) = UCONV_E(K) |
| 565 |
C VGS(J,K) = VCONV_S(K) |
| 566 |
C VGN(J,K) = VCONV_N(K) |
| 567 |
END DO |
| 568 |
|
| 569 |
|
| 570 |
|
| 571 |
c update the cloud-base mass flux: |
| 572 |
|
| 573 |
CBMFG1(J) = CBMF |
| 574 |
PLCLG1(J) = PLCL |
| 575 |
|
| 576 |
c fill-in the MITPHYGR2 common: |
| 577 |
|
| 578 |
DO K = 1, NL |
| 579 |
MUP1(J,K) = MUP(K) |
| 580 |
MDOWN1(J,K) = MDOWN(K) |
| 581 |
MP1(J,K) = MP0(K) |
| 582 |
ENT1(J,K) = ENT(K) |
| 583 |
DET1(J,K) = DET(K) |
| 584 |
ENDDO |
| 585 |
INBG1(J) = INB |
| 586 |
ICBG1(J) = ICB |
| 587 |
|
| 588 |
c ------------------------------------------------------------------------- |
| 589 |
c 3. Cloudiness parameterization |
| 590 |
c |
| 591 |
c (the subroutine predicts the cloud fraction and cloud water associated |
| 592 |
c with cumulus convection, and the precipitation and the cloud water |
| 593 |
c associated with large-scale super-saturation) |
| 594 |
c |
| 595 |
c CLDF: cloud fraction [0-1] (will be used in radiation) |
| 596 |
c |
| 597 |
c CLDQ: in-cloud mixing ratio of condensed (liquid+solid) water [kg/kg] |
| 598 |
c resulting from cumulus convection + large-scale condensation |
| 599 |
c (will be used to compute cloud optical properties). |
| 600 |
c ------------------------------------------------------------------------- |
| 601 |
|
| 602 |
|
| 603 |
C modif omp: does not compute the cloudfraction if CLEARSKY is .TRUE. |
| 604 |
|
| 605 |
IF (CLEARSKY) THEN |
| 606 |
DO K = 1,NL |
| 607 |
QCONDC(K) = 0. |
| 608 |
END DO |
| 609 |
END IF |
| 610 |
|
| 611 |
CALL CLOUDS_SUB_LS_40(NLEV,RCONV,RS,TCONV,PPLAY,PAPRS |
| 612 |
: ,DELT,QCONDC |
| 613 |
: ,CLDF,CLDQ,PRADJ,FTADJ,FRADJ) |
| 614 |
|
| 615 |
|
| 616 |
C modif omp. 12-30-01 (Some people do not take enough vacations...) |
| 617 |
C --> average the cloud over between the radiation calls. |
| 618 |
C --> average for the cloud fraction and total water. |
| 619 |
|
| 620 |
DO K = 1,NLEV |
| 621 |
CLDF_MEAN(J,K) = CLDF_MEAN(J,K) + CLDF(K) * DELT |
| 622 |
CLDQ_MEAN(J,K) = CLDQ_MEAN(J,K) + CLDF(K) * CLDQ(K) * DELT |
| 623 |
END DO |
| 624 |
|
| 625 |
c T and Q tendencies: |
| 626 |
|
| 627 |
DO K = 1, NLEV |
| 628 |
TT_LSC(J,K) = FTADJ(K) |
| 629 |
QT_LSC(J,K) = FRADJ(K) |
| 630 |
ENDDO |
| 631 |
|
| 632 |
c ----------------------------------------------------------------------- |
| 633 |
c 4. Cloud optical properties |
| 634 |
c |
| 635 |
c (interface between clouds and radiation). |
| 636 |
c |
| 637 |
c The following variables will be used in input of the radiation code: |
| 638 |
c |
| 639 |
c CLDF: cloud fraction [0-1] |
| 640 |
c CLDEMI: cloud longwave emissivity [0-1] |
| 641 |
c CLDTAU: cloud visible optical thickness |
| 642 |
c ----------------------------------------------------------------------- |
| 643 |
|
| 644 |
DO K = 1, NLEV |
| 645 |
PPLAY1(K) = PPLAY(K)*100. |
| 646 |
PAPRS1(K) = PAPRS(K)*100. |
| 647 |
ENDDO |
| 648 |
PAPRS1(NLEV+1) = PAPRS(NLEV+1)*100. |
| 649 |
|
| 650 |
IF ( RADINT ) |
| 651 |
: CALL OPTICAL(NLEV,TCONV,PAPRS1,PPLAY1,CLDF,CLDQ |
| 652 |
: ,CLDEMI,CLDTAU,CLDFICE,CLDT,CLDWP) |
| 653 |
|
| 654 |
c ----------------------------------------------------------------------- |
| 655 |
c 5. Radiation |
| 656 |
c |
| 657 |
c Important: to save some CPU time, radiation is not computed at every |
| 658 |
c timestep but every RADFREQ seconds. |
| 659 |
c ----------------------------------------------------------------------- |
| 660 |
|
| 661 |
IF ( RADINT ) THEN ! interactive radiation |
| 662 |
|
| 663 |
c If desired, call radiation scheme |
| 664 |
|
| 665 |
C IF ( CURTIME .LT. (1.5*DELT) .OR. |
| 666 |
C : RADTIME .GT. (RADFREQ-10.0) .OR. Firstcall ) THEN |
| 667 |
|
| 668 |
IF ( (RADTIME .GT. (RADFREQ-10.0)) .OR. Firstcall ) THEN |
| 669 |
|
| 670 |
IF (J .EQ. NGP) then |
| 671 |
write(*,*) 'Radiation call: time = ',CURTIME |
| 672 |
write(*,*) 'julien =', julien |
| 673 |
RADTIME=0.0 |
| 674 |
END IF |
| 675 |
|
| 676 |
|
| 677 |
C modif omp: mean cloud cover |
| 678 |
|
| 679 |
DO K = 1,NLEV |
| 680 |
CLDF_M(K) = CLDF_MEAN(J,K) / RADFREQ |
| 681 |
IF (CLDF_M(K) .GT. 0) THEN |
| 682 |
CLDQ_M(K) = CLDQ_MEAN(J,K) / RADFREQ / CLDF_M(K) |
| 683 |
ELSE |
| 684 |
CLDQ_M(K) = 0. |
| 685 |
END IF |
| 686 |
|
| 687 |
CLDF_MEAN(J,K) = 0.0 |
| 688 |
CLDQ_MEAN(J,K) = 0.0 |
| 689 |
END DO |
| 690 |
CALL OPTICAL(NLEV,TCONV,PAPRS1,PPLAY1,CLDF_M,CLDQ_M |
| 691 |
& ,CLDEMI_M,CLDTAU_M,CLDFICE_M,CLDT_M,CLDWP_M) |
| 692 |
|
| 693 |
c IF (DARAD.eq.'n'.or.DARAD.eq.'N') THEN |
| 694 |
c TIMEU=TIME0+MOD(CURTIME,86400) |
| 695 |
c ELSE |
| 696 |
c TIMEU=TIME0 |
| 697 |
c ENDIF |
| 698 |
|
| 699 |
C fixed time for radiation |
| 700 |
|
| 701 |
C TIMEU = 24. * 3600. * 31. * 28. * 20. ! current time in seconds |
| 702 |
TIMEU = rad_start_time ! current time in seconds |
| 703 |
|
| 704 |
C adjust time from starting time |
| 705 |
IF (RAD_TIME) THEN |
| 706 |
TIMEU = TIMEU +CURTIME ! current time in seconds |
| 707 |
END IF |
| 708 |
|
| 709 |
IF (RAD_LAT) THEN |
| 710 |
RLAT = LAT_G(J) ! latitude |
| 711 |
ELSE |
| 712 |
RLAT = 0. ! latitude (valid only for equatorial model...) |
| 713 |
END IF |
| 714 |
IF (RAD_LON) THEN |
| 715 |
RLON = LON_G(J) |
| 716 |
ELSE |
| 717 |
RLON = 0. ! longitude (here we want a uniform insolation) |
| 718 |
END IF |
| 719 |
c -- Earth-Sun distance: |
| 720 |
|
| 721 |
julien = INT(TIMEU)/(3600*24) + 1 |
| 722 |
|
| 723 |
c omp: Classic fortran bug - single precision real passed to a subroutine |
| 724 |
C requiring a doouble precision. |
| 725 |
C do the conversion independently of the precision of REAL... |
| 726 |
rjulien = julien |
| 727 |
|
| 728 |
C CALL ORBITE(FLOAT(julien),zlongi,dist) |
| 729 |
C CALL ORBITE(DFLOAT(julien),zlongi,dist) |
| 730 |
CALL ORBITE(rjulien,zlongi,dist) |
| 731 |
|
| 732 |
c -- solar zenith angle and surface albedo: |
| 733 |
|
| 734 |
C IF(DARAD.EQ.'y'.OR.DARAD.EQ.'Y')THEN ! no diurnal cycle |
| 735 |
IF( DARAD )THEN ! no diurnal cycle |
| 736 |
CALL ANGLE(zlongi,RLAT,fract,rmu0) |
| 737 |
C omp: same thing here. |
| 738 |
C CALL ALBOC(FLOAT(julien),RLAT,albsol) |
| 739 |
C CALL ALBOC(DFLOAT(julien),RLAT,albsol) |
| 740 |
CALL ALBOC(rjulien,RLAT,albsol) |
| 741 |
albsol = ALB ! only if you prefer to set it to a constant value |
| 742 |
ELSE |
| 743 |
gmtime = MOD(TIMEU,3600.*24.)/86400. - RADFREQ/2./86400. |
| 744 |
CALL ZENANG(zlongi,gmtime,RADFREQ,RLAT,RLON,rmu0,fract) |
| 745 |
CALL ALBOC_CD(rmu0,albsol) |
| 746 |
ENDIF |
| 747 |
|
| 748 |
c -- call the radiation subroutine |
| 749 |
|
| 750 |
TSRAD = TS(J) |
| 751 |
|
| 752 |
C modif omp: mean cloud cover |
| 753 |
|
| 754 |
DO K = 1, NLEV |
| 755 |
C instantaneous value of the cloud cover |
| 756 |
RVCONV(K) = RCONV(K) - CLDQ(K)*CLDF(K) |
| 757 |
C Mean value |
| 758 |
C RVCONV(K) = RCONV(K) - CLDQ_M(K)*CLDF_M(K) |
| 759 |
|
| 760 |
RVCONV(K) = MAX( MIN(RVCONV(K),RCONV(K)) , 0.0 ) |
| 761 |
O3RAD(K) = O3(K) * 1.0E-6 |
| 762 |
|
| 763 |
IF (CLEARSKY) THEN |
| 764 |
CLDFRAD(K) = 0. |
| 765 |
CLDEMIRAD(K) = 0. |
| 766 |
CLDTAURAD(K) = 0. |
| 767 |
ELSE |
| 768 |
C instantaneous value of the cloud cover |
| 769 |
CLDFRAD(K) = CLDF(K) |
| 770 |
CLDEMIRAD(K) = CLDEMI(K) |
| 771 |
CLDTAURAD(K) = CLDTAU(K) |
| 772 |
C mean value for the cloud cover |
| 773 |
C CLDFRAD(K) = CLDF_M(K) |
| 774 |
C CLDEMIRAD(K) = CLDEMI_M(K) |
| 775 |
C CLDTAURAD(K) = CLDTAU_M(K) |
| 776 |
ENDIF |
| 777 |
|
| 778 |
ENDDO |
| 779 |
|
| 780 |
CALL RADLWSW |
| 781 |
i (dist, rmu0, fract, CO2, SCON, |
| 782 |
i PAPRS1, PPLAY1,TSRAD,albsol,TCONV,RVCONV,O3RAD, |
| 783 |
i CLDFRAD, CLDEMIRAD, CLDTAURAD, |
| 784 |
o heat,heat0,cool,cool0,radsol,albpla, |
| 785 |
o topsw,toplw,solsw,sollw, |
| 786 |
o topsw0,toplw0,solsw0,sollw0,inso) |
| 787 |
|
| 788 |
DO K = 1, NLEV |
| 789 |
|
| 790 |
c if thermo constants used in radiation are different than those |
| 791 |
c used in the calling program, renormalize the cooling rates: |
| 792 |
|
| 793 |
normalis = 1004.709/(CPD*(1.-RCONV(K))+CPV*RCONV(K)) |
| 794 |
: *G/9.80665 |
| 795 |
cool(K) = cool(K)*normalis |
| 796 |
heat(K) = heat(K)*normalis |
| 797 |
cool0(K) = cool0(K)*normalis |
| 798 |
heat0(K) = heat0(K)*normalis |
| 799 |
|
| 800 |
ENDDO |
| 801 |
|
| 802 |
TLW(J) = toplw |
| 803 |
TSW(J) = topsw |
| 804 |
TLW0(J) = toplw0 |
| 805 |
TSW0(J) = topsw0 |
| 806 |
SLW(J) = sollw |
| 807 |
SSW(J) = solsw |
| 808 |
SLW0(J) = sollw0 |
| 809 |
SSW0(J) = solsw0 |
| 810 |
SINS(J) = inso |
| 811 |
|
| 812 |
c T and Q tendencies associated with radiation: |
| 813 |
|
| 814 |
IF ( CRFINT ) THEN |
| 815 |
|
| 816 |
DO K = 1, NLEV |
| 817 |
TT_RLW(J,K) = -cool(K)/(24*3600.) |
| 818 |
TT_RSW(J,K) = heat(K)/(24*3600.) |
| 819 |
TT_RLW0(J,K) = -cool0(K)/(24*3600.) |
| 820 |
TT_RSW0(J,K) = heat0(K)/(24*3600.) |
| 821 |
|
| 822 |
RC1(J,K) = heat(K)-cool(K) |
| 823 |
CRLW1(J,K) = cool0(K)-cool(K) |
| 824 |
CRSW1(J,K) = heat(K)-heat0(K) |
| 825 |
ENDDO |
| 826 |
|
| 827 |
ELSE ! CRFINT=false |
| 828 |
|
| 829 |
c caution: cloudy rad fluxes at the surface and at the top of the atmosphere |
| 830 |
c are not meaningfull if the CRF is specified! (but clear-sky fluxes are) |
| 831 |
|
| 832 |
DO K = 1, NLEV |
| 833 |
c clear-sky rates are those actually computed: |
| 834 |
TT_RLW0(J,K) = -cool0(K)/(24*3600.) |
| 835 |
TT_RSW0(J,K) = heat0(K)/(24*3600.) |
| 836 |
c cloudy rates deduced from "clear-sky+specified CRF": |
| 837 |
c be careful: the names of the following variables may not correspond |
| 838 |
c to their actual definition here (misleading notations in that case): |
| 839 |
RC1(J,K) = heat0(K)-cool0(K) + RCRF0(K) ! net rad heating rate |
| 840 |
TT_RLW(J,K) = (heat0(K)-cool0(K)+RCRF0(K))/(24*3600.)!net rad heating rate here |
| 841 |
TT_RSW(J,K) = 0.0 ! anything (will not be used) |
| 842 |
CRLW1(J,K) = RCRF0(K) ! here: net CRF (notation misleading) |
| 843 |
CRSW1(J,K) = heat0(K)-cool0(K) ! here: net clear-sky heating rate (misleading) |
| 844 |
ENDDO |
| 845 |
|
| 846 |
|
| 847 |
ENDIF ! CRFINT |
| 848 |
|
| 849 |
ENDIF ! RADFREQ |
| 850 |
|
| 851 |
ELSE ! RADINT = FALSE |
| 852 |
|
| 853 |
DO K = 1,NLEV |
| 854 |
|
| 855 |
TT_RLW(J,K) = RC1(J,K)/86400. |
| 856 |
TT_RSW(J,K) = 0.0 |
| 857 |
TT_RLW0(J,K) = 0.0 |
| 858 |
TT_RSW0(J,K) = 0.0 |
| 859 |
CRLW1(J,K) = 0.0 |
| 860 |
CRSW1(J,K) = 0.0 |
| 861 |
END DO |
| 862 |
TLW(J) = 0.0 |
| 863 |
TSW(J) = 0.0 |
| 864 |
TLW0(J) = 0.0 |
| 865 |
TSW0(J) = 0.0 |
| 866 |
SLW(J) = 0.0 |
| 867 |
SSW(J) = 0.0 |
| 868 |
SLW0(J) = 0.0 |
| 869 |
SSW0(J) = 0.0 |
| 870 |
SINS(J) = 0.0 |
| 871 |
|
| 872 |
ENDIF ! RADINT |
| 873 |
|
| 874 |
|
| 875 |
c ----------------------------------------------------------------------- |
| 876 |
c 6. Surface fluxes: |
| 877 |
c ----------------------------------------------------------------------- |
| 878 |
|
| 879 |
C omp: mixed-layer |
| 880 |
IF (PBL_MIX) THEN |
| 881 |
|
| 882 |
U_PBLW = 0.0 |
| 883 |
U_PBLE = 0.0 |
| 884 |
V_PBLS = 0.0 |
| 885 |
V_PBLN = 0.0 |
| 886 |
|
| 887 |
PBL_INV = 1./float(K_PBL) |
| 888 |
DO k = 1,K_PBL |
| 889 |
U_PBLW = U_PBLW + UGW(J,k) |
| 890 |
U_PBLE = U_PBLE + UGE(J,k) |
| 891 |
V_PBLS = V_PBLS + VGS(J,k) |
| 892 |
V_PBLN = V_PBLN + VGN(J,k) |
| 893 |
END DO |
| 894 |
U_PBLW = U_PBLW * PBL_INV |
| 895 |
U_PBLE = U_PBLE * PBL_INV |
| 896 |
V_PBLS = V_PBLS * PBL_INV |
| 897 |
V_PBLN = V_PBLN * PBL_INV |
| 898 |
|
| 899 |
C OMP - This version is incorrect. |
| 900 |
C The vertical velocities must also be modified in order for |
| 901 |
C advection (which is computed after the call to MITPHYS) |
| 902 |
C to be computed properly. |
| 903 |
C DO k = 1,K_PBL |
| 904 |
C UGW(J,K) = U_PBLW |
| 905 |
C UGE(J,K) = U_PBLE |
| 906 |
C VGS(J,K) = V_PBLS |
| 907 |
C VGN(J,K) = V_PBLN |
| 908 |
C END DO |
| 909 |
C Alternative: include the change in the _PBL tendencies: |
| 910 |
|
| 911 |
|
| 912 |
DO k = 1, K_PBL |
| 913 |
UT_PBL_E (J,K) = 0.5 * DELTI_MIX * (U_PBLE- UGE(J,K)) |
| 914 |
UT_PBL_W (J,K) = 0.5 * DELTI_MIX * (U_PBLW- UGW(J,K)) |
| 915 |
VT_PBL_S (J,K) = 0.5 * DELTI_MIX * (V_PBLS- VGS(J,K)) |
| 916 |
VT_PBL_N (J,K) = 0.5 * DELTI_MIX * (V_PBLN- VGN(J,K)) |
| 917 |
END DO |
| 918 |
|
| 919 |
|
| 920 |
END IF |
| 921 |
|
| 922 |
|
| 923 |
C omp: older version |
| 924 |
C Vertical diffusion in the PBL instead of mixing. It does not work as well |
| 925 |
C when the PBL is deep, and still allows significant vertical shear. |
| 926 |
|
| 927 |
c$$$ VISC_PBL = .1 *((paprs1(1) - paprs1(2))) **2 / DELT |
| 928 |
c$$$C GAMMA = RD_AIR / CP_AIR |
| 929 |
c$$$ |
| 930 |
c$$$ DPI = 1./(pplay1(K-1) - pplay1(K) ) |
| 931 |
c$$$ |
| 932 |
c$$$ DO K = 2, K_PBL |
| 933 |
c$$$ UE_FLUX(K) = VISC_PBL * |
| 934 |
c$$$ & ( UGE(J,K-1) - UGE(J,K)) * DPI |
| 935 |
c$$$ UW_FLUX(K) = VISC_PBL * |
| 936 |
c$$$ & ( UGW(J,K-1) - UGW(J,K)) * DPI |
| 937 |
c$$$ VS_FLUX(K) = VISC_PBL * |
| 938 |
c$$$ & ( VGS(J,K-1) - VGS(J,K)) * DPI |
| 939 |
c$$$ VN_FLUX(K) = VISC_PBL * |
| 940 |
c$$$ & ( VGN(J,K-1) - VGN(J,K)) * DPI |
| 941 |
c$$$ |
| 942 |
c$$$C QT_FLUX(K) = VISC_PBL * |
| 943 |
c$$$C & ( QG1(J,K-1) - QG1(J,K)) * DPI |
| 944 |
c$$$C QT_FLUX(K) = 0. |
| 945 |
c$$$C H_FLUX(K) = 0. |
| 946 |
c$$$C H_FLUX(K) = VISC_PBL * |
| 947 |
c$$$C & ( TG1(I,J,K-1) * ( P_full(K-1) ** GAMMA) |
| 948 |
c$$$C & - TG1(I,J,K) * ( P_full(K) ** GAMMA)) |
| 949 |
c$$$C & * (P_half(K) ** (- GAMMA)) * DPI |
| 950 |
c$$$ END DO |
| 951 |
c$$$C QT_FLUX(I_PBL) = 0. |
| 952 |
c$$$C H_FLUX(I_PBL) = 0. |
| 953 |
c$$$ UW_FLUX(1) = 0. |
| 954 |
c$$$ UE_FLUX(1) = 0. |
| 955 |
c$$$ VS_FLUX(1) = 0. |
| 956 |
c$$$ VN_FLUX(1) = 0. |
| 957 |
c$$$ UW_FLUX(K_PBL+1) = 0. |
| 958 |
c$$$ UE_FLUX(K_PBL+1) = 0. |
| 959 |
c$$$ VS_FLUX(K_PBL+1) = 0. |
| 960 |
c$$$ VN_FLUX(K_PBL+1) = 0. |
| 961 |
c$$$ DO K = 1,K_PBL |
| 962 |
c$$$ DPI = 1./(paprs1(K) - paprs1(K+1) ) |
| 963 |
c$$$C QT_PBL(I,J,K) = (QT_FLUX(K) - QT_FLUX(K+1)) |
| 964 |
c$$$C & * DPI |
| 965 |
c$$$C TT_PBL(I,J,K) = (H_FLUX(K) - H_FLUX(K+1)) |
| 966 |
c$$$C & * DPI |
| 967 |
c$$$ UT_PBL_E(J,K) = 0.5 * (UE_FLUX(K) - UE_FLUX(K+1)) |
| 968 |
c$$$ & * DPI |
| 969 |
c$$$ UT_PBL_W(J,K) = 0.5 * (UW_FLUX(K) - UW_FLUX(K+1)) |
| 970 |
c$$$ & * DPI |
| 971 |
c$$$ VT_PBL_S(J,K) = 0.5 *(VS_FLUX(K) - VS_FLUX(K+1)) |
| 972 |
c$$$ & * DPI |
| 973 |
c$$$ VT_PBL_N(J,K) = 0.5 *(VN_FLUX(K) - VN_FLUX(K+1)) |
| 974 |
c$$$ & * DPI |
| 975 |
c$$$ END DO |
| 976 |
c$$$ |
| 977 |
c$$$ |
| 978 |
c$$$ |
| 979 |
|
| 980 |
c -- surface wind: |
| 981 |
c US0, VS0: mean surface wind |
| 982 |
c WS0: minimum surface wind |
| 983 |
c WD: gustiness factor due to convective downdrafts |
| 984 |
|
| 985 |
|
| 986 |
c omp: two formulation for the contribution of the zonal wind: |
| 987 |
C |
| 988 |
|
| 989 |
C 1) US2 and VS2 are the square of the velocity at the center of the |
| 990 |
C cell |
| 991 |
|
| 992 |
US2= 0.25 * (US0+LOG_CORRECT *UGW(J,1) |
| 993 |
& + US0+LOG_CORRECT *UGE(J,1)) * |
| 994 |
& (US0+LOG_CORRECT *UGW(J,1) + |
| 995 |
& US0+LOG_CORRECT *UGE(J,1)) |
| 996 |
|
| 997 |
VS2= 0.25 *( VS0+LOG_CORRECT *VGS(J,1) |
| 998 |
& + VS0+LOG_CORRECT *VGN(J,1)) * |
| 999 |
& (VS0+LOG_CORRECT *VGS(J,1) |
| 1000 |
& + VS0+LOG_CORRECT *VGN(J,1)) |
| 1001 |
|
| 1002 |
|
| 1003 |
C or 2) US2 and VS2 are the average of the square of the velocity |
| 1004 |
|
| 1005 |
C US2= 0.5 *( (US0+LOG_CORRECT *UGW(J,1))* |
| 1006 |
C & (US0+LOG_CORRECT *UGW(J,1)) |
| 1007 |
C & + (US0+LOG_CORRECT *UGE(J,1)) |
| 1008 |
C & * (US0+LOG_CORRECT *UGE(J,1))) |
| 1009 |
|
| 1010 |
|
| 1011 |
C VS2= 0.5 *( (VS0+LOG_CORRECT *VGS(J,1))* |
| 1012 |
C & (VS0+LOG_CORRECT *VGS(J,1)) |
| 1013 |
C & + (VS0+LOG_CORRECT *VGN(J,1)) |
| 1014 |
C & * (VS0+LOG_CORRECT *VGN(J,1))) |
| 1015 |
|
| 1016 |
VSURF=SQRT(WS0*WS0+WD*WD+VS2+US2) |
| 1017 |
VSPRIME=VSURF-SQRT(WS0*WS0+VS2+US2) |
| 1018 |
|
| 1019 |
|
| 1020 |
c -- turbulent fluxes (NB: PG in mb): |
| 1021 |
|
| 1022 |
TSA=TCONV(1)*(PG/PPLAY(1))**(RD/CPD) |
| 1023 |
ROWS=PG/(RD*TSA*(1.+RCONV(1)*(EPSI-1.))) |
| 1024 |
DPH1I=1.0/(PAPRS(1)-PAPRS(2)) |
| 1025 |
|
| 1026 |
FTSURF = G * ROWS * CD |
| 1027 |
: * (VSURF*(TS(J)-TSA) - VSPRIME*TPRIME)*DPH1I |
| 1028 |
|
| 1029 |
FRSURF = G * ROWS * CD * EFRAC * LANDMUL(J) |
| 1030 |
: * (VSURF*(RSS-RCONV(1)) - VSPRIME*QPRIME)*DPH1I |
| 1031 |
|
| 1032 |
c -- sensible and latent heat fluxes at the surface (W/m2): |
| 1033 |
|
| 1034 |
TC = TS(J)-273.15 |
| 1035 |
ALV = LV0-CPVMCL*TC |
| 1036 |
|
| 1037 |
SH = (CPD*(1.-RCONV(1))+CPV*RCONV(1)) * 100.0*ROWS * CD |
| 1038 |
: * (VSURF*(TS(J)-TSA) - VSPRIME*TPRIME) |
| 1039 |
|
| 1040 |
LH = ALV * 100.0*ROWS * CD * EFRAC * LANDMUL(J) |
| 1041 |
: * (VSURF*(RSS-RCONV(1)) - VSPRIME*QPRIME) |
| 1042 |
|
| 1043 |
c -- Surface flux of momentum: |
| 1044 |
c here we impose the mean flow and so we apply |
| 1045 |
c the drag only to the wind perturbation: |
| 1046 |
|
| 1047 |
IF (LINEAR_FRICTION) THEN |
| 1048 |
|
| 1049 |
VSURF = LIN_FR_VEL |
| 1050 |
END IF |
| 1051 |
FUSURF_E = -G*ROWS*CDM*VSURF*DPH1I*UGE(J,1) * LOG_CORRECT |
| 1052 |
FUSURF_W = -G*ROWS*CDM*VSURF*DPH1I*UGW(J,1) * LOG_CORRECT |
| 1053 |
|
| 1054 |
C FUSURF_E = -G*ROWS*CDM*VSURF*DPH1I* 0.5 |
| 1055 |
C & * (UGE(J,1) + UGW(J,1)) * LOG_CORRECT |
| 1056 |
C FUSURF_W = FUSURF_E |
| 1057 |
|
| 1058 |
|
| 1059 |
FVSURF_S = -G*ROWS*CDM*VSURF*DPH1I*VGS(J,1) * LOG_CORRECT |
| 1060 |
FVSURF_N = -G*ROWS*CDM*VSURF*DPH1I*VGN(J,1) * LOG_CORRECT |
| 1061 |
|
| 1062 |
C FVSURF_S = -G*ROWS*CDM*VSURF*DPH1I*0.5 |
| 1063 |
C & * (VGS(J,1) + VGN(J,1)) * LOG_CORRECT |
| 1064 |
C FVSURF_N = FVSURF_S |
| 1065 |
|
| 1066 |
|
| 1067 |
c -- relaxation of the zonally averaged zonal wind |
| 1068 |
c to a preset value: |
| 1069 |
|
| 1070 |
c DO K = 1, NLEV |
| 1071 |
c UT_PBLi(J,K) = (UMEAN-UZON(K)) / TAUM |
| 1072 |
c VT_PBL(J,K) = 0. |
| 1073 |
c ENDDO |
| 1074 |
|
| 1075 |
c -- T, Q, U, V tendencies at first atmospheric level: |
| 1076 |
|
| 1077 |
TT_PBL(J,1) = FTSURF |
| 1078 |
QT_PBL(J,1) = FRSURF |
| 1079 |
UT_PBL_E(J,1) = UT_PBL_E(J,1) + 0.5 * FUSURF_E |
| 1080 |
UT_PBL_W(J,1) = UT_PBL_W(J,1) + 0.5 * FUSURF_W |
| 1081 |
VT_PBL_S(J,1) = VT_PBL_S(J,1) + 0.5 * FVSURF_S |
| 1082 |
VT_PBL_N(J,1) = VT_PBL_N(J,1) + 0.5 * FVSURF_N |
| 1083 |
|
| 1084 |
|
| 1085 |
c ----------------------------------------------------------------------- |
| 1086 |
c 7. Fill-in COMMONs : |
| 1087 |
c ----------------------------------------------------------------------- |
| 1088 |
|
| 1089 |
c MITFLUXES: |
| 1090 |
|
| 1091 |
PRECNV(J) = PRECIP ! convective precipitation (mm/day) |
| 1092 |
PRECLS(J) = PRADJ ! large-scale precipitation (mm/day) |
| 1093 |
LHF(J) = LH |
| 1094 |
SHF(J) = SH |
| 1095 |
|
| 1096 |
USTR(J) = 0.5 * (FUSURF_E + FUSURF_W) |
| 1097 |
& * (PAPRS(1)-PAPRS(2)) / G * 100.0 |
| 1098 |
VSTR(J) = 0.5 * (FVSURF_S + FVSURF_N) |
| 1099 |
& * (PAPRS(1)-PAPRS(2)) / G * 100.0 |
| 1100 |
|
| 1101 |
|
| 1102 |
c MITPHYGR3: |
| 1103 |
|
| 1104 |
CLDT1(J) = CLDT ! total cloud cover |
| 1105 |
CLDWP1(J) = CLDWP ! cloud water path |
| 1106 |
TS1(J) = TS(J) ! surface temperature |
| 1107 |
CC (acz) check for STL |
| 1108 |
CC IF( J .EQ. NGP ) THEN |
| 1109 |
CC write(6,*) 'TS1 set to TS: TS1=',TS1(J) |
| 1110 |
CC END IF |
| 1111 |
PRW1(J) = 0.0 |
| 1112 |
DO K = 1, NLEV |
| 1113 |
AUX = (RCONV(K)-CLDF(K)*CLDQ(K))/RS(K) |
| 1114 |
RH1(J,K) = MIN(MAX(AUX,0.0),1.0) |
| 1115 |
CLDF1(J,K) = CLDF(K) |
| 1116 |
CLDQ1(J,K) = CLDQ(K) |
| 1117 |
CLDQC1(J,K) = QCONDC(K) |
| 1118 |
PRW1(J) = PRW1(J)+RCONV(K)*(PAPRS1(K)-PAPRS1(K+1))/G |
| 1119 |
ENDDO |
| 1120 |
|
| 1121 |
GZ(1)=0.0 |
| 1122 |
CPN(1)=CPD*(1.-RCONV(1))+RCONV(1)*CPV |
| 1123 |
LV(1)=LV0-CPVMCL*(TCONV(1)-273.15) |
| 1124 |
SX(1)=CPN(I)*TCONV(1)+GZ(1) |
| 1125 |
HX(1)=SX(1)+LV(1)*RCONV(1) |
| 1126 |
DO K = 2, NLEV |
| 1127 |
CPN(K)=CPD*(1.-RCONV(K))+RCONV(K)*CPV |
| 1128 |
TVX=TCONV(K)*(1.+RCONV(K)*EPSI-RCONV(K)) |
| 1129 |
TVY=TCONV(K-1)*(1.+RCONV(K-1)*EPSI-RCONV(K-1)) |
| 1130 |
GZ(K)=GZ(K-1)+0.5*RD*(TVX+TVY) |
| 1131 |
: *(PPLAY(K-1)-PPLAY(K))/PAPRS(K) |
| 1132 |
LV(K)=LV0-CPVMCL*(TCONV(K)-273.15) |
| 1133 |
SX(K)=TCONV(K)*CPN(K)+GZ(K) |
| 1134 |
HX(K)=SX(K)+LV(K)*RCONV(K) |
| 1135 |
ENDDO |
| 1136 |
DO K = 1, NLEV |
| 1137 |
SXG1(J,K)=SX(K)/1000. |
| 1138 |
HXG1(J,K)=HX(K)/1000. |
| 1139 |
ENDDO |
| 1140 |
|
| 1141 |
CC (acz) interactive SST and STL: computed at each |
| 1142 |
CC gridpoint (then masked to determine TS) |
| 1143 |
CC with different QFLUX (defined in do_inphys) |
| 1144 |
CC |
| 1145 |
IF (.NOT. PRESC_SST) THEN |
| 1146 |
SST_TEND = (-SH - LH + SSW(J) - SLW(J) + QFLUX(J)) |
| 1147 |
: / HEATCAPOCEA |
| 1148 |
SST1(J) = SST1(J) + SST_TEND * DELT |
| 1149 |
END IF |
| 1150 |
|
| 1151 |
IF (.NOT. PRESC_STL) THEN |
| 1152 |
STL_TEND = (-SH - LH + SSW(J) - SLW(J) + QFLUX(J)) |
| 1153 |
: / HEATCAPLAND |
| 1154 |
STL1(J) = STL1(J) + STL_TEND * DELT |
| 1155 |
END IF |
| 1156 |
|
| 1157 |
CC |
| 1158 |
CC (nprive+acz) Interactive soil moisture |
| 1159 |
CC |
| 1160 |
|
| 1161 |
c convert latent heat flux to evaporation and |
| 1162 |
c evaporation in kg m-2 s-1 into evaporation in cm/s |
| 1163 |
EBKT = LH / (10. * ALV) |
| 1164 |
|
| 1165 |
c convert precipitation in mm/day in cm/s |
| 1166 |
PBKT = (PRECIP + PRADJ)/(10*24*3600) |
| 1167 |
|
| 1168 |
IF (.NOT. PRESC_BKT) THEN |
| 1169 |
BUCKET(J) = BUCKET(J) + (PBKT-EBKT) * DELT |
| 1170 |
IF (BUCKET(J) .LT. 0) then |
| 1171 |
BUCKET(J) = 0 |
| 1172 |
ELSEIF (BUCKET(J) .GT. BKT0) then |
| 1173 |
BUCKET(J) = BKT0 |
| 1174 |
ENDIF |
| 1175 |
|
| 1176 |
ELSE |
| 1177 |
|
| 1178 |
BUCKET(J) = BKT0 |
| 1179 |
|
| 1180 |
ENDIF |
| 1181 |
|
| 1182 |
|
| 1183 |
C ================================================================ |
| 1184 |
|
| 1185 |
9999 CONTINUE ! loop over atmospheric columns ends up here |
| 1186 |
|
| 1187 |
C ================================================================ |
| 1188 |
|
| 1189 |
C DUMMY = 0 |
| 1190 |
C DO k = 1, K_PBL |
| 1191 |
C DUMMY = DUMMY + UT_PBL_W(100,k) + UT_PBL_E(100,k) |
| 1192 |
C END DO |
| 1193 |
|
| 1194 |
C write(*,*) 'MITPHYS_DRIVER:' |
| 1195 |
C write(*,*) 'UT_PBL(100,1:K_PBL):', DUMMY |
| 1196 |
C write(*,*) 'USTR(100) = ',USTR(100) |
| 1197 |
C DUMMY = 0 |
| 1198 |
C DO k = 1, K_PBL |
| 1199 |
C DUMMY = DUMMY + VT_PBL_S(100,k) + VT_PBL_N(100,k) |
| 1200 |
C END DO |
| 1201 |
|
| 1202 |
C write(*,*) 'VT_PBL(100,1):', DUMMY |
| 1203 |
C write(*,*) 'VSTR(100) = ', VSTR(100) |
| 1204 |
C write(*,*) ROWS |
| 1205 |
C |
| 1206 |
C write(*,*) 'PHY_DRIVER end, TG1(65,1):', TG1(65,1) |
| 1207 |
|
| 1208 |
C write(*,*) 'PHY_DRIVER end, CBMFG1(85):', CBMFG1(85) |
| 1209 |
|
| 1210 |
FirstCall = .FALSE. |
| 1211 |
|
| 1212 |
RETURN |
| 1213 |
END |
| 1214 |
|
| 1215 |
|
| 1216 |
C $Id: phy_driver.F,v 1.2 2000/06/26 23:45:42 cnh Exp $ |
| 1217 |
|
| 1218 |
|
| 1219 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph