#include "CPP_OPTIONS.h" SUBROUTINE MPDRIVER (CURTIME, DELT) IMPLICIT NONE C-- C-- SUBROUTINE PDRIVER (CURTIME) C-- C-- Purpose: stand-alone driver for MIT physical parametrization routines C-- Input : CURTIME : current time (in seconds) C-- grid-point model fields in common block: PHYGR1 C-- forcing fields in common blocks : LSMASK, FORFIX, FORCIN C-- Output : Diagnosed convective variables in common block: MITPHYGR2 C-- Other diagnosed variables in common block: MITPHYGR3 C-- Physical param. tendencies in common block: MITPHYTEN C-- Surface and upper boundary fluxes in common block: MITFLUXES C-- C C============================================================================= C C ******************************************************************* C * * C * MIT Physics package * C * * C * Version 1.0 * C * (July 2000) * C * * C ******************************************************************* C C C Created: 07/2000 Sandrine Bony & Kerry Emanuel C Modification by Olivier Pauluis (07-2001 --- now) C C -- Important: C ------------- C C * this version (1.0) of the physics package is intended to be used over C a tropical aqua-planet (no sea ice, no distinction between land and C ocean surfaces, no specific cloudiness associated with mid-latitude C weather systems, etc.) C ... but nothing precludes to complete it later on ! C C * for a non-equatorial version of the model, don't forget to specify C the latitude (RLAT) seen by radiation. C C * if you change the vertical resolution, change it accordingly C in the file ../inc/dimensions.h (this is for the radiation code) C and also in convect43b.F (NLEV=...) C (not very satisfying...will have to be improved) C C C -- Among important differences with Franco Molteni's physics: C ------------------------------------------------------------- C C - opposite convention used for numbering the vertical levels C (here, the 1st atmospheric level is the nearest to the surface) C C - QG1 (specific humidity) is in kg/kg C C - all the physics package is called atmospheric column by C atmospheric column (note that the radiation code would C easily accept to work with 3D fields) C C C -- If you are looking for some documentation about the parameterizations C used in this package: C -------------------------------------------------------------------------- C C Bony S and K A Emanuel, 2001: A parameterization of the cloudiness C associated with cumulus convection; Evaluation using TOGA COARE C data. J. Atmos. Sci., accepted. C C Emanuel K A, 1991: A scheme for representing cumulus convection in C large-scale models. J. Atmos. Sci., 48, 2313-2335. C C Emanuel K A and M Zivkovic-Rothman, 1999: Development and evaluation C of a convection scheme for use in climate models. J. Atmos. Sci., C 56, 1766-1782 (see also http://www-paoc.mit.edu/~emanuel/home.html) C C Fouquart Y and B Bonnel, 1980: Computation of solar heating of the C Earth's atmosphere: a new parameterization. Beitr. Phys. Atmos., C 53, 35-62. C C Morcrette J-J, 1991: Radiation and cloud radiative properties in the C European Centre for Medium-Range Weather Forecasts forecasting C system. J. Geophys. Res., 96, 9121-9132. C C-- Please, kindly report any problem or suggestion to the authors: C C bony@wind.mit.edu, emanuel@texmex.mit.edu C C============================================================================= C Resolution parameters C #include "atparam.h" #include "atparam1.h" #include "MITPHYS_PARAMS.h" C INTEGER NLON, NLAT, NLEV, NGP PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT ) C C Constants + functions of sigma and latitude C #include "Lev_def.h" C C Model variables, tendencies and fluxes on gaussian grid C #include "com_mitphysvar.h" C C Diagnostics from the convection scheme: C #include "com_convect.h" C C Assign thermodynamical constants: C #include "thermo.h" C C Surface forcing fields (time-inv. or functions of seasonal cycle) C #include "com_forcing1.h" c Misc: INTEGER I,J,K _RL CURTIME _RL DELT, ES, AUX, RSS, RS(NLEV), QMIN PARAMETER (QMIN=1.e-12) ! minimum humidity when negative C PARAMETER (DELT=600.D0) ! timestep in seconds (temporary: deltaT better) C (not very satisfying...will have to be improved) _RL TS(NGP) _RL SX(NLEV),HX(NLEV),LV(NLEV),GZ(NLEV),TVY,TVX,CPN(NLEV) c Convection: INTEGER IFLAG, NTRA, NL PARAMETER (NTRA=1) ! number of tracors in convection _RL CBMF, PRECIP, WD, TPRIME, QPRIME _RL TCONV(NLEV), RCONV(NLEV), RVCONV(NLEV) _RL TMEMO(NLEV), RMEMO(NLEV) _RL UCONV_E(NLEV), UCONV_W(NLEV) _RL VCONV_S(NLEV), VCONV_N(NLEV) _RL UCONV(NLEV), VCONV(NLEV) _RL TRACONV(NLEV,NTRA) _RL QCONDC(NLEV), FTCONV(NLEV), FRCONV(NLEV) _RL FUCONV(NLEV), FVCONV(NLEV), FTRACONV(NLEV,NTRA) C modif omp _RL PLCL c Soil moisture: bucket model _RL LANDMUL(NGP), FACTOR, EBKT, PBKT c Cloudiness and large-scale condensation: _RL PRADJ, CLDT, CLDWP _RL CLDF(NLEV), CLDQ(NLEV), FTADJ(NLEV), FRADJ(NLEV) _RL CLDEMI(NLEV), CLDTAU(NLEV), CLDFICE(NLEV) _RL CLDFRAD(NLEV), CLDEMIRAD(NLEV), CLDTAURAD(NLEV) _RL CLDT_M, CLDWP_M _RL CLDF_M(NLEV), CLDQ_M(NLEV) _RL CLDEMI_M(NLEV), CLDTAU_M(NLEV), CLDFICE_M(NLEV) c Radiation: c modif omp: now included in MITPHYS_PARAMS.h C _RL RADFREQ ! frequency of radiation calls C PARAMETER (RADFREQ=6*3600.) ! here: every 6 hours _RL SCON, CO2, ALB! solar cst, CO2 conc, sfc albedo PARAMETER (SCON=1370.0, CO2=330.0) !units: W/m2, ppm c PARAMETER (ALB=0.07) !units: fraction PARAMETER (ALB=0.10) !units: fraction C REPLACED by a logical DARAD in MITPHYS_PARAMS.h C CHARACTER*1 DARAD ! diurnally averaged radiation? C PARAMETER (DARAD='y') _RL RC0(NGP,NLEV) ! specified radiative heating rate _RL RCRF0(NLEV) ! specified net cld radiative forcing C modif omp: these parameters are defined in the MITPHYS namelist C LOGICAL RADINT, CRFINT, CLEARSKY C PARAMETER (RADINT=.TRUE.) ! TRUE=interactive radiation C PARAMETER (CRFINT=.TRUE.)! TRUE=interactive cld rad forcing C PARAMETER (CLEARSKY=.FALSE.)! TRUE=ignore cld-rad interactions C PARAMETER (RADINT=.FALSE.) ! TRUE=interactive radiation C PARAMETER (CRFINT=.FALSE.)! TRUE=interactive cld rad forcing C PARAMETER (CLEARSKY=.TRUE.)! TRUE=ignore cld-rad interactions INTEGER julien _RL rjulien _RL zlongi, dist, rmu0, fract, TSRAD, TIMEU, gmtime _RL radsol, albsol, albpla, RLON, RLAT, normalis _RL O3(NLEV), O3RAD(NLEV) _RL topsw, toplw, solsw, sollw _RL topsw0, toplw0, solsw0, sollw0, inso _RL heat(NLEV), heat0(NLEV), cool(NLEV), cool0(NLEV) _RL RADTIME SAVE RADTIME LOGICAL FirstCall DATA FirstCall /.TRUE./ SAVE FirstCall c Surface fluxes: C modif omp: these parameters are defined in the MITPHYS namelist C _RL CD ! drag coeff and sfc evap potential C PARAMETER (CD=1.2D-3) C _RL US0,VS0,WS0 ! mean and min sfc winds C PARAMETER (US0= 0.0,VS0=0.0,WS0=4.0) _RL EFRAC PARAMETER (EFRAC = 1.0) _RL US2,VS2,VSURF,VSPRIME,TC,ALV _RL TSA,PG,ROWS,DPH1I,FTSURF,FRSURF,SH,LH _RL SST_TEND, STL_TEND _RL DUMMY C modif omp: these parameters are defined in the MITPHYS namelist C OMP 07-18-01 C LOG_CORRECT: correction factor for wind distribution in the PBL. C Defined as the ratio of the 10m wind to the wind at the C lowest model level C LOG_CORRECT = 1: wind uniform in the lowest level C LOG_CORRECT = 2/3: logarithmic profile between the surface to 100m, with a C surfcae roughness ~ 10cm. C _RL LOG_CORRECT C PARAMETER ( LOG_CORRECT = 1.0 ) C PARAMETER ( LOG_CORRECT = 0.6666666666 ) C OMP: LOW LEVEL MIXING: C modif omp: these parameters are defined in the MITPHYS namelist C INTEGER K_PBL C PARAMETER (K_PBL = 8 ) C LOGICAL PBL_DIFF C PARAMETER (PBL_DIFF = .FALSE.) _RL U_PBLW, U_PBLE, V_PBLS, V_PBLN, PBL_INV _RL UW_FLUX(nlev), UE_FLUX(nlev), VN_FLUX(nlev),VS_FLUX(nlev) _RL DPI, VISC_PBL, DELTI_ADJ, DELTI_MIX _RL WIND_ADJ _RL FUSURF_W, FUSURF_E, FVSURF_S,FVSURF_N c Momentum drag: c _RL TAUM ! timescale for zonal wind relaxation c PARAMETER (TAUM=3.*86400.) ! 3 days _RL UMEAN, VMEAN ! specified mean barotropic flow PARAMETER (UMEAN=0.0, VMEAN=0.) ! for drag only C modif omp: these parameters are defined in the MITPHYS namelist C _RL CDM ! drag coeff for sfc flux of momentum C PARAMETER (CDM=1.2d-3) c _RL UZON(NLEV) C Atmospheric pressure levels used by the physics package: C PPLAY (mb): pressure surfaces where T, Q are computed C PAPRS (mb): pressure surfaces at the interface between layers C PPLAY1 and PAPRS1: same as PPLAY and PAPRS but in Pa C (should be improved later...use pSurfs instead) _RL DP, PPLAY(NLEV), PAPRS(NLEV+1) _RL PPLAY1(NLEV), PAPRS1(NLEV+1) c 5levels DATA PPLAY / 950.0, 775.0, 500.0, 250.0, 75.0 / c 5levels DATA PAPRS / 1000.0, 900.0, 650.0, 350.0, 150.0, 0.0 / c 5levels DATA O3 / 0.0, 0.0, 0.0, 0.0, 0.0 / DATA O3 / 0.048, 0.049, 0.051, 0.052, 0.053, 0.054, 0.056, : 0.056, 0.057, 0.058, 0.059, 0.059, 0.060, 0.060, : 0.061, 0.062, 0.063, 0.064, 0.066, 0.068, 0.069, : 0.072, 0.073, 0.075, 0.077, 0.081, 0.085, 0.090, : 0.096, 0.106, 0.117, 0.135, 0.153, 0.175, 0.202, : 0.243, 0.405, 1.000, 3.700, 10.124 / c mean CRF profile derived for SST=27C, U*=2m/s, U0=-5m/s, 360x40, 180E: DATA RCRF0 / 0.152041972, 0.183760509, 0.167694584, 0.16433309, : 0.16215764, 0.161808997, 0.163289174, 0.166284353, : 0.170786828, 0.176764682, 0.1845745, 0.19467622, : 0.201630697, 0.209879801, 0.223531589, 0.245427981, : 0.270941138, 0.301004857, 0.324049205, 0.340905219, : 0.360249668, 0.389321774, 0.412838638, 0.449643821, : 0.483938843, 0.525295079, 0.666303575, 0.935836792, : 1.25989616, 1.19505048, 1.40494013, 1.28129566, : -0.570507288, -0.167021349, -0.00548044546, 0.129761115, : -0.0882121176, -0.0795640424, -0.134504586, 0.0340689607 / C ========================================================================= C C Loop over atmospheric columns starts here: C C ========================================================================= CC (acz) check for SST CC write(6,*) 'beginning of mit_phydriver: SST1=',SST1(1) IF (FirstCall) THEN DO J = 1,NGP DO K = 1 , NLEV CLDF_MEAN(J,K) = 0.0 CLDQ_MEAN(J,K) = 0.0 END DO END DO END IF C note: using 1.0/DELT potentially unstable because of the C Adams-Bashford time-stepping. CAN be put to 1 if the physics is C taken out of the Adams-Bashford... IF (CONV_ADJ) THEN DELTI_ADJ = 1.0/ DELT WIND_ADJ = 0.6 ELSE DELTI_ADJ = 1.0/ DELT WIND_ADJ = 0.0 END IF DELTI_MIX = 0.6 / DELT RADTIME = RADTIME + DELT c -- pressure levels: DP = 1000./NLEV DO K = 1, NLEV PPLAY(K) = 1000. - (K-1)*DP IF(K.GT.1) PAPRS(K)=0.5*(PPLAY(K)+PPLAY(K-1)) ENDDO PAPRS(1) = 2.*PPLAY(1)-PAPRS(2) PAPRS(NLEV+1) = 2.*PAPRS(NLEV)-PAPRS(NLEV-1) PAPRS(NLEV+1) = MAX(0.1,PAPRS(NLEV+1)) c -- zonally-averaged zonal wind at each level: c DO K = 1, NLEV c UZON(K) = 0. c DO J = 1, NGP c UZON(K) = UZON(K) + UG1(J,K)/FLOAT(NGP) c ENDDO c ENDDO c -- loop over horizontal gridpoints: C Newtonian cooling rate IF (.NOT. RADINT ) THEN CALL NEWT_COOL(TG1,RC0) END IF DO 9999 J = 1, NGP c ----------------------------------------------------------------------- c 1. Preliminaries: c ----------------------------------------------------------------------- c -- check that humidities are not negative: DO K = 1, NLEV if (QG1(J,K).lt.0.0 ) then CC(acz) write(*,*) ' J, K, neg humidity: ',J, K,QG1(J,K) QG1(J,K) = MAX( QG1(J,K), QMIN ) endif ENDDO c -- initialize tendencies (MITPHYTEN common): DO K = 1, NLEV TT_CNV(J,K) = 0.0 QT_CNV(J,K) = 0.0 UT_CNV(J,K) = 0.0 VT_CNV(J,K) = 0.0 TT_LSC(J,K) = 0.0 QT_LSC(J,K) = 0.0 TT_PBL(J,K) = 0.0 QT_PBL(J,K) = 0.0 UT_PBL_E(J,K) = 0.0 UT_PBL_W(J,K) = 0.0 VT_PBL_S(J,K) = 0.0 VT_PBL_N(J,K) = 0.0 ENDDO PRECNV(J) = 0. PRECLS(J) = 0. LHF(J) = 0. SHF(J) = 0. c -- surface pressure (mb): PG = 0.01 * PSG1(J) c -- surface temperature (K): TS(J) =SST1(J)+FMASK1(J)*(STL1(J)-SST1(J)) CC (acz) check for SST CC write(6,*) 'start of phy_driver: J, FMASK1(J)=',J, FMASK1(J) CC if( J .EQ. NGP) THEN CC write(6,*) 'start of phy_driver: TS, SST1, STL1=',TS(NGP), SST1(NGP),STL1(NGP) CC end if c -- compute the land multiplicator factor for evaporation c (acz) adapted from (nprive) - Delworth and Manabe (1988) if (BUCKET(J) .GE. 0.75*BKT0) then FACTOR = 1 else FACTOR = BUCKET(J)/(0.75*BKT0) endif LANDMUL(J) = 1 + FMASK1(J) * (FACTOR-1) c -- compute saturation specific humidity: TC=TS(J)-273.15 ES=6.112*EXP(17.67*TC/(243.5+TC)) RSS=0.98*EPS*ES/(PG-(1.-EPS)*ES) DO K = 1, NLEV TC=TG1(J,K)-273.15 IF(TC.GE.0.0)THEN ES=6.112*EXP(17.67*TC/(243.5+TC)) ELSE ES=EXP(23.33086-6111.72784/TG1(J,K)+0.15215*LOG(TG1(J,K))) ENDIF ES=MIN(ES,PPLAY(K)) RS(K)=EPS*ES/(PPLAY(K)-(1.-EPS)*ES) ENDDO c -- eventually, prescribe the radiative cooling rate to a cst value: C omp: Prescribed cooling has been replaced by a newtonian cooling IF ( .NOT. RADINT ) THEN DO K = 1, NLEV C RC1(J,K) = - RC0(J,K)! net rad heating (units = K/day) RC1(J,K) = RC0(J,K) * 24.d0 * 3600.d0 ! net rad heating (units = K/day) ENDDO ENDIF c ----------------------------------------------------------------------- c 2. Convection subroutine c c CBMF : cloud-base mass flux (its value must be "remembered" by the c calling program between calls to the convection subroutine). c c QCONDC: in-cloud mixing ratio of condensed water produced by cumulus c convection [kg/kg], output of the convection scheme used in c the cloud parameterization. c c WD, TPRIME, QPRIME: convective downdraft velocity scale [m/s], c temperature perturbation scale [K] and specific humidity c perturbation scale [kg/kg] that will be used in surface c flux parameterizations. c ----------------------------------------------------------------------- c cloud-base mass flux: CBMF = CBMFG1(J) PLCL = 0 c NL: the maximum number of levels to which convection can penetrate, plus 1 c (NL MUST be less than or equal to NLEV-1) NL = 1 DO K = 1, NLEV IF(PPLAY(K).GE.78.0) NL = MAX(NL,K) ENDDO c T and Q profiles seen by the convection (and eventually altered by c the subroutine is dry convective adjustment occurs): DO K = 1, NLEV TCONV(K) = TG1(J,K) ! TCONV in K RCONV(K) = QG1(J,K) ! RCONV in kg/kg C OMP: Current implementation of the momentum transport consider the mean C velocity at the center of the cell. Future implementatioin may C separate between UGW, VGW, VGS, VGN. this would require modifying C convect43 to handle more than 2 velocity. UCONV(K) = 0.5 * (UGW(J,K) + UGE(J,K)) VCONV(K) = 0.5 * (VGS(J,K) + VGN(J,K)) UCONV_W(K) = UGW(J,K) UCONV_E(K) = UGE(J,K) VCONV_S(K) = VGS(J,K) VCONV_N(K) = VGN(J,K) DO I = 1, NTRA TRACONV(K,I) = 0.0 ! no tracors here ENDDO ENDDO c "official" version 4.3b of CONVECT plus diagnostics required for cld parameterization: C write(*,*) '--- convect: I2 = ', J, 'CBMF =', CBMF CALL CONVECT43C(TCONV, RCONV, RS, UCONV, VCONV, : TRACONV, : PPLAY, PAPRS, NLEV, NL, NTRA, DELT, : IFLAG, FTCONV, FRCONV, FUCONV, FVCONV, : FTRACONV, PRECIP, WD, TPRIME, QPRIME, QCONDC, : CBMF, PLCL ) C write(*,*) 'after convect, CBMF =', CBMF C if necessary, write to error file IF (IFLAG.EQ.4) THEN WRITE(*,120) CURTIME 120 FORMAT(5X,'CFL condition violated at TIME ',e15.7) WRITE(*,*) 'LAT: ', LAT_G(J), 'LON: ', LON_G(J) ENDIF c T and Q tendencies associated with convection: DO K = 1, NL TT_CNV(J,K) = FTCONV(K) ! K/s QT_CNV(J,K) = FRCONV(K) ! kg/kg/s UT_CNV(J,K) = FUCONV(K) ! m/s^2 VT_CNV(J,K) = FVCONV(K) ! m/s^2 c c note that we do not update the wind, here c (should be done in principle if the dry adjustment c is called and if the wind-tendencies due to the c convection are to be taken into account). c modif omp: pass the convective adjustment as a time tendecy TT_ADJ(J,K) = DELTI_ADJ * (TCONV(K) - TG1(J,K)) QT_ADJ(J,K) = DELTI_ADJ * (RCONV(K) - QG1(J,K)) UT_ADJ(J,K) = WIND_ADJ * DELTI_ADJ * (UCONV(K) & - 0.5 * (UCONV_W(K) + UCONV_E(K))) VT_ADJ(J,K) = WIND_ADJ * DELTI_ADJ * (VCONV(K) & - 0.5 * (VCONV_S(K) + VCONV_N(K))) TG1(J,K) = TCONV(K) QG1(J,K) = RCONV(K) ! QG1 in kg/kg C UGW(J,K) = UCONV_W(K) C UGE(J,K) = UCONV_E(K) C VGS(J,K) = VCONV_S(K) C VGN(J,K) = VCONV_N(K) END DO c update the cloud-base mass flux: CBMFG1(J) = CBMF PLCLG1(J) = PLCL c fill-in the MITPHYGR2 common: DO K = 1, NL MUP1(J,K) = MUP(K) MDOWN1(J,K) = MDOWN(K) MP1(J,K) = MP0(K) ENT1(J,K) = ENT(K) DET1(J,K) = DET(K) ENDDO INBG1(J) = INB ICBG1(J) = ICB c ------------------------------------------------------------------------- c 3. Cloudiness parameterization c c (the subroutine predicts the cloud fraction and cloud water associated c with cumulus convection, and the precipitation and the cloud water c associated with large-scale super-saturation) c c CLDF: cloud fraction [0-1] (will be used in radiation) c c CLDQ: in-cloud mixing ratio of condensed (liquid+solid) water [kg/kg] c resulting from cumulus convection + large-scale condensation c (will be used to compute cloud optical properties). c ------------------------------------------------------------------------- C modif omp: does not compute the cloudfraction if CLEARSKY is .TRUE. IF (CLEARSKY) THEN DO K = 1,NL QCONDC(K) = 0. END DO END IF CALL CLOUDS_SUB_LS_40(NLEV,RCONV,RS,TCONV,PPLAY,PAPRS : ,DELT,QCONDC : ,CLDF,CLDQ,PRADJ,FTADJ,FRADJ) C modif omp. 12-30-01 (Some people do not take enough vacations...) C --> average the cloud over between the radiation calls. C --> average for the cloud fraction and total water. DO K = 1,NLEV CLDF_MEAN(J,K) = CLDF_MEAN(J,K) + CLDF(K) * DELT CLDQ_MEAN(J,K) = CLDQ_MEAN(J,K) + CLDF(K) * CLDQ(K) * DELT END DO c T and Q tendencies: DO K = 1, NLEV TT_LSC(J,K) = FTADJ(K) QT_LSC(J,K) = FRADJ(K) ENDDO c ----------------------------------------------------------------------- c 4. Cloud optical properties c c (interface between clouds and radiation). c c The following variables will be used in input of the radiation code: c c CLDF: cloud fraction [0-1] c CLDEMI: cloud longwave emissivity [0-1] c CLDTAU: cloud visible optical thickness c ----------------------------------------------------------------------- DO K = 1, NLEV PPLAY1(K) = PPLAY(K)*100. PAPRS1(K) = PAPRS(K)*100. ENDDO PAPRS1(NLEV+1) = PAPRS(NLEV+1)*100. IF ( RADINT ) : CALL OPTICAL(NLEV,TCONV,PAPRS1,PPLAY1,CLDF,CLDQ : ,CLDEMI,CLDTAU,CLDFICE,CLDT,CLDWP) c ----------------------------------------------------------------------- c 5. Radiation c c Important: to save some CPU time, radiation is not computed at every c timestep but every RADFREQ seconds. c ----------------------------------------------------------------------- IF ( RADINT ) THEN ! interactive radiation c If desired, call radiation scheme C IF ( CURTIME .LT. (1.5*DELT) .OR. C : RADTIME .GT. (RADFREQ-10.0) .OR. Firstcall ) THEN IF ( (RADTIME .GT. (RADFREQ-10.0)) .OR. Firstcall ) THEN IF (J .EQ. NGP) then write(*,*) 'Radiation call: time = ',CURTIME write(*,*) 'julien =', julien RADTIME=0.0 END IF C modif omp: mean cloud cover DO K = 1,NLEV CLDF_M(K) = CLDF_MEAN(J,K) / RADFREQ IF (CLDF_M(K) .GT. 0) THEN CLDQ_M(K) = CLDQ_MEAN(J,K) / RADFREQ / CLDF_M(K) ELSE CLDQ_M(K) = 0. END IF CLDF_MEAN(J,K) = 0.0 CLDQ_MEAN(J,K) = 0.0 END DO CALL OPTICAL(NLEV,TCONV,PAPRS1,PPLAY1,CLDF_M,CLDQ_M & ,CLDEMI_M,CLDTAU_M,CLDFICE_M,CLDT_M,CLDWP_M) c IF (DARAD.eq.'n'.or.DARAD.eq.'N') THEN c TIMEU=TIME0+MOD(CURTIME,86400) c ELSE c TIMEU=TIME0 c ENDIF C fixed time for radiation C TIMEU = 24. * 3600. * 31. * 28. * 20. ! current time in seconds TIMEU = rad_start_time ! current time in seconds C adjust time from starting time IF (RAD_TIME) THEN TIMEU = TIMEU +CURTIME ! current time in seconds END IF IF (RAD_LAT) THEN RLAT = LAT_G(J) ! latitude ELSE RLAT = 0. ! latitude (valid only for equatorial model...) END IF IF (RAD_LON) THEN RLON = LON_G(J) ELSE RLON = 0. ! longitude (here we want a uniform insolation) END IF c -- Earth-Sun distance: julien = INT(TIMEU)/(3600*24) + 1 c omp: Classic fortran bug - single precision real passed to a subroutine C requiring a doouble precision. C do the conversion independently of the precision of REAL... rjulien = julien C CALL ORBITE(FLOAT(julien),zlongi,dist) C CALL ORBITE(DFLOAT(julien),zlongi,dist) CALL ORBITE(rjulien,zlongi,dist) c -- solar zenith angle and surface albedo: C IF(DARAD.EQ.'y'.OR.DARAD.EQ.'Y')THEN ! no diurnal cycle IF( DARAD )THEN ! no diurnal cycle CALL ANGLE(zlongi,RLAT,fract,rmu0) C omp: same thing here. C CALL ALBOC(FLOAT(julien),RLAT,albsol) C CALL ALBOC(DFLOAT(julien),RLAT,albsol) CALL ALBOC(rjulien,RLAT,albsol) albsol = ALB ! only if you prefer to set it to a constant value ELSE gmtime = MOD(TIMEU,3600.*24.)/86400. - RADFREQ/2./86400. CALL ZENANG(zlongi,gmtime,RADFREQ,RLAT,RLON,rmu0,fract) CALL ALBOC_CD(rmu0,albsol) ENDIF c -- call the radiation subroutine TSRAD = TS(J) C modif omp: mean cloud cover DO K = 1, NLEV C instantaneous value of the cloud cover RVCONV(K) = RCONV(K) - CLDQ(K)*CLDF(K) C Mean value C RVCONV(K) = RCONV(K) - CLDQ_M(K)*CLDF_M(K) RVCONV(K) = MAX( MIN(RVCONV(K),RCONV(K)) , 0.0 ) O3RAD(K) = O3(K) * 1.0E-6 IF (CLEARSKY) THEN CLDFRAD(K) = 0. CLDEMIRAD(K) = 0. CLDTAURAD(K) = 0. ELSE C instantaneous value of the cloud cover CLDFRAD(K) = CLDF(K) CLDEMIRAD(K) = CLDEMI(K) CLDTAURAD(K) = CLDTAU(K) C mean value for the cloud cover C CLDFRAD(K) = CLDF_M(K) C CLDEMIRAD(K) = CLDEMI_M(K) C CLDTAURAD(K) = CLDTAU_M(K) ENDIF ENDDO CALL RADLWSW i (dist, rmu0, fract, CO2, SCON, i PAPRS1, PPLAY1,TSRAD,albsol,TCONV,RVCONV,O3RAD, i CLDFRAD, CLDEMIRAD, CLDTAURAD, o heat,heat0,cool,cool0,radsol,albpla, o topsw,toplw,solsw,sollw, o topsw0,toplw0,solsw0,sollw0,inso) DO K = 1, NLEV c if thermo constants used in radiation are different than those c used in the calling program, renormalize the cooling rates: normalis = 1004.709/(CPD*(1.-RCONV(K))+CPV*RCONV(K)) : *G/9.80665 cool(K) = cool(K)*normalis heat(K) = heat(K)*normalis cool0(K) = cool0(K)*normalis heat0(K) = heat0(K)*normalis ENDDO TLW(J) = toplw TSW(J) = topsw TLW0(J) = toplw0 TSW0(J) = topsw0 SLW(J) = sollw SSW(J) = solsw SLW0(J) = sollw0 SSW0(J) = solsw0 SINS(J) = inso c T and Q tendencies associated with radiation: IF ( CRFINT ) THEN DO K = 1, NLEV TT_RLW(J,K) = -cool(K)/(24*3600.) TT_RSW(J,K) = heat(K)/(24*3600.) TT_RLW0(J,K) = -cool0(K)/(24*3600.) TT_RSW0(J,K) = heat0(K)/(24*3600.) RC1(J,K) = heat(K)-cool(K) CRLW1(J,K) = cool0(K)-cool(K) CRSW1(J,K) = heat(K)-heat0(K) ENDDO ELSE ! CRFINT=false c caution: cloudy rad fluxes at the surface and at the top of the atmosphere c are not meaningfull if the CRF is specified! (but clear-sky fluxes are) DO K = 1, NLEV c clear-sky rates are those actually computed: TT_RLW0(J,K) = -cool0(K)/(24*3600.) TT_RSW0(J,K) = heat0(K)/(24*3600.) c cloudy rates deduced from "clear-sky+specified CRF": c be careful: the names of the following variables may not correspond c to their actual definition here (misleading notations in that case): RC1(J,K) = heat0(K)-cool0(K) + RCRF0(K) ! net rad heating rate TT_RLW(J,K) = (heat0(K)-cool0(K)+RCRF0(K))/(24*3600.)!net rad heating rate here TT_RSW(J,K) = 0.0 ! anything (will not be used) CRLW1(J,K) = RCRF0(K) ! here: net CRF (notation misleading) CRSW1(J,K) = heat0(K)-cool0(K) ! here: net clear-sky heating rate (misleading) ENDDO ENDIF ! CRFINT ENDIF ! RADFREQ ELSE ! RADINT = FALSE DO K = 1,NLEV TT_RLW(J,K) = RC1(J,K)/86400. TT_RSW(J,K) = 0.0 TT_RLW0(J,K) = 0.0 TT_RSW0(J,K) = 0.0 CRLW1(J,K) = 0.0 CRSW1(J,K) = 0.0 END DO TLW(J) = 0.0 TSW(J) = 0.0 TLW0(J) = 0.0 TSW0(J) = 0.0 SLW(J) = 0.0 SSW(J) = 0.0 SLW0(J) = 0.0 SSW0(J) = 0.0 SINS(J) = 0.0 ENDIF ! RADINT c ----------------------------------------------------------------------- c 6. Surface fluxes: c ----------------------------------------------------------------------- C omp: mixed-layer IF (PBL_MIX) THEN U_PBLW = 0.0 U_PBLE = 0.0 V_PBLS = 0.0 V_PBLN = 0.0 PBL_INV = 1./float(K_PBL) DO k = 1,K_PBL U_PBLW = U_PBLW + UGW(J,k) U_PBLE = U_PBLE + UGE(J,k) V_PBLS = V_PBLS + VGS(J,k) V_PBLN = V_PBLN + VGN(J,k) END DO U_PBLW = U_PBLW * PBL_INV U_PBLE = U_PBLE * PBL_INV V_PBLS = V_PBLS * PBL_INV V_PBLN = V_PBLN * PBL_INV C OMP - This version is incorrect. C The vertical velocities must also be modified in order for C advection (which is computed after the call to MITPHYS) C to be computed properly. C DO k = 1,K_PBL C UGW(J,K) = U_PBLW C UGE(J,K) = U_PBLE C VGS(J,K) = V_PBLS C VGN(J,K) = V_PBLN C END DO C Alternative: include the change in the _PBL tendencies: DO k = 1, K_PBL UT_PBL_E (J,K) = 0.5 * DELTI_MIX * (U_PBLE- UGE(J,K)) UT_PBL_W (J,K) = 0.5 * DELTI_MIX * (U_PBLW- UGW(J,K)) VT_PBL_S (J,K) = 0.5 * DELTI_MIX * (V_PBLS- VGS(J,K)) VT_PBL_N (J,K) = 0.5 * DELTI_MIX * (V_PBLN- VGN(J,K)) END DO END IF C omp: older version C Vertical diffusion in the PBL instead of mixing. It does not work as well C when the PBL is deep, and still allows significant vertical shear. c$$$ VISC_PBL = .1 *((paprs1(1) - paprs1(2))) **2 / DELT c$$$C GAMMA = RD_AIR / CP_AIR c$$$ c$$$ DPI = 1./(pplay1(K-1) - pplay1(K) ) c$$$ c$$$ DO K = 2, K_PBL c$$$ UE_FLUX(K) = VISC_PBL * c$$$ & ( UGE(J,K-1) - UGE(J,K)) * DPI c$$$ UW_FLUX(K) = VISC_PBL * c$$$ & ( UGW(J,K-1) - UGW(J,K)) * DPI c$$$ VS_FLUX(K) = VISC_PBL * c$$$ & ( VGS(J,K-1) - VGS(J,K)) * DPI c$$$ VN_FLUX(K) = VISC_PBL * c$$$ & ( VGN(J,K-1) - VGN(J,K)) * DPI c$$$ c$$$C QT_FLUX(K) = VISC_PBL * c$$$C & ( QG1(J,K-1) - QG1(J,K)) * DPI c$$$C QT_FLUX(K) = 0. c$$$C H_FLUX(K) = 0. c$$$C H_FLUX(K) = VISC_PBL * c$$$C & ( TG1(I,J,K-1) * ( P_full(K-1) ** GAMMA) c$$$C & - TG1(I,J,K) * ( P_full(K) ** GAMMA)) c$$$C & * (P_half(K) ** (- GAMMA)) * DPI c$$$ END DO c$$$C QT_FLUX(I_PBL) = 0. c$$$C H_FLUX(I_PBL) = 0. c$$$ UW_FLUX(1) = 0. c$$$ UE_FLUX(1) = 0. c$$$ VS_FLUX(1) = 0. c$$$ VN_FLUX(1) = 0. c$$$ UW_FLUX(K_PBL+1) = 0. c$$$ UE_FLUX(K_PBL+1) = 0. c$$$ VS_FLUX(K_PBL+1) = 0. c$$$ VN_FLUX(K_PBL+1) = 0. c$$$ DO K = 1,K_PBL c$$$ DPI = 1./(paprs1(K) - paprs1(K+1) ) c$$$C QT_PBL(I,J,K) = (QT_FLUX(K) - QT_FLUX(K+1)) c$$$C & * DPI c$$$C TT_PBL(I,J,K) = (H_FLUX(K) - H_FLUX(K+1)) c$$$C & * DPI c$$$ UT_PBL_E(J,K) = 0.5 * (UE_FLUX(K) - UE_FLUX(K+1)) c$$$ & * DPI c$$$ UT_PBL_W(J,K) = 0.5 * (UW_FLUX(K) - UW_FLUX(K+1)) c$$$ & * DPI c$$$ VT_PBL_S(J,K) = 0.5 *(VS_FLUX(K) - VS_FLUX(K+1)) c$$$ & * DPI c$$$ VT_PBL_N(J,K) = 0.5 *(VN_FLUX(K) - VN_FLUX(K+1)) c$$$ & * DPI c$$$ END DO c$$$ c$$$ c$$$ c -- surface wind: c US0, VS0: mean surface wind c WS0: minimum surface wind c WD: gustiness factor due to convective downdrafts c omp: two formulation for the contribution of the zonal wind: C C 1) US2 and VS2 are the square of the velocity at the center of the C cell US2= 0.25 * (US0+LOG_CORRECT *UGW(J,1) & + US0+LOG_CORRECT *UGE(J,1)) * & (US0+LOG_CORRECT *UGW(J,1) + & US0+LOG_CORRECT *UGE(J,1)) VS2= 0.25 *( VS0+LOG_CORRECT *VGS(J,1) & + VS0+LOG_CORRECT *VGN(J,1)) * & (VS0+LOG_CORRECT *VGS(J,1) & + VS0+LOG_CORRECT *VGN(J,1)) C or 2) US2 and VS2 are the average of the square of the velocity C US2= 0.5 *( (US0+LOG_CORRECT *UGW(J,1))* C & (US0+LOG_CORRECT *UGW(J,1)) C & + (US0+LOG_CORRECT *UGE(J,1)) C & * (US0+LOG_CORRECT *UGE(J,1))) C VS2= 0.5 *( (VS0+LOG_CORRECT *VGS(J,1))* C & (VS0+LOG_CORRECT *VGS(J,1)) C & + (VS0+LOG_CORRECT *VGN(J,1)) C & * (VS0+LOG_CORRECT *VGN(J,1))) VSURF=SQRT(WS0*WS0+WD*WD+VS2+US2) VSPRIME=VSURF-SQRT(WS0*WS0+VS2+US2) c -- turbulent fluxes (NB: PG in mb): TSA=TCONV(1)*(PG/PPLAY(1))**(RD/CPD) ROWS=PG/(RD*TSA*(1.+RCONV(1)*(EPSI-1.))) DPH1I=1.0/(PAPRS(1)-PAPRS(2)) FTSURF = G * ROWS * CD : * (VSURF*(TS(J)-TSA) - VSPRIME*TPRIME)*DPH1I FRSURF = G * ROWS * CD * EFRAC * LANDMUL(J) : * (VSURF*(RSS-RCONV(1)) - VSPRIME*QPRIME)*DPH1I c -- sensible and latent heat fluxes at the surface (W/m2): TC = TS(J)-273.15 ALV = LV0-CPVMCL*TC SH = (CPD*(1.-RCONV(1))+CPV*RCONV(1)) * 100.0*ROWS * CD : * (VSURF*(TS(J)-TSA) - VSPRIME*TPRIME) LH = ALV * 100.0*ROWS * CD * EFRAC * LANDMUL(J) : * (VSURF*(RSS-RCONV(1)) - VSPRIME*QPRIME) c -- Surface flux of momentum: c here we impose the mean flow and so we apply c the drag only to the wind perturbation: IF (LINEAR_FRICTION) THEN VSURF = LIN_FR_VEL END IF FUSURF_E = -G*ROWS*CDM*VSURF*DPH1I*UGE(J,1) * LOG_CORRECT FUSURF_W = -G*ROWS*CDM*VSURF*DPH1I*UGW(J,1) * LOG_CORRECT C FUSURF_E = -G*ROWS*CDM*VSURF*DPH1I* 0.5 C & * (UGE(J,1) + UGW(J,1)) * LOG_CORRECT C FUSURF_W = FUSURF_E FVSURF_S = -G*ROWS*CDM*VSURF*DPH1I*VGS(J,1) * LOG_CORRECT FVSURF_N = -G*ROWS*CDM*VSURF*DPH1I*VGN(J,1) * LOG_CORRECT C FVSURF_S = -G*ROWS*CDM*VSURF*DPH1I*0.5 C & * (VGS(J,1) + VGN(J,1)) * LOG_CORRECT C FVSURF_N = FVSURF_S c -- relaxation of the zonally averaged zonal wind c to a preset value: c DO K = 1, NLEV c UT_PBLi(J,K) = (UMEAN-UZON(K)) / TAUM c VT_PBL(J,K) = 0. c ENDDO c -- T, Q, U, V tendencies at first atmospheric level: TT_PBL(J,1) = FTSURF QT_PBL(J,1) = FRSURF UT_PBL_E(J,1) = UT_PBL_E(J,1) + 0.5 * FUSURF_E UT_PBL_W(J,1) = UT_PBL_W(J,1) + 0.5 * FUSURF_W VT_PBL_S(J,1) = VT_PBL_S(J,1) + 0.5 * FVSURF_S VT_PBL_N(J,1) = VT_PBL_N(J,1) + 0.5 * FVSURF_N c ----------------------------------------------------------------------- c 7. Fill-in COMMONs : c ----------------------------------------------------------------------- c MITFLUXES: PRECNV(J) = PRECIP ! convective precipitation (mm/day) PRECLS(J) = PRADJ ! large-scale precipitation (mm/day) LHF(J) = LH SHF(J) = SH USTR(J) = 0.5 * (FUSURF_E + FUSURF_W) & * (PAPRS(1)-PAPRS(2)) / G * 100.0 VSTR(J) = 0.5 * (FVSURF_S + FVSURF_N) & * (PAPRS(1)-PAPRS(2)) / G * 100.0 c MITPHYGR3: CLDT1(J) = CLDT ! total cloud cover CLDWP1(J) = CLDWP ! cloud water path TS1(J) = TS(J) ! surface temperature CC (acz) check for STL CC IF( J .EQ. NGP ) THEN CC write(6,*) 'TS1 set to TS: TS1=',TS1(J) CC END IF PRW1(J) = 0.0 DO K = 1, NLEV AUX = (RCONV(K)-CLDF(K)*CLDQ(K))/RS(K) RH1(J,K) = MIN(MAX(AUX,0.0),1.0) CLDF1(J,K) = CLDF(K) CLDQ1(J,K) = CLDQ(K) CLDQC1(J,K) = QCONDC(K) PRW1(J) = PRW1(J)+RCONV(K)*(PAPRS1(K)-PAPRS1(K+1))/G ENDDO GZ(1)=0.0 CPN(1)=CPD*(1.-RCONV(1))+RCONV(1)*CPV LV(1)=LV0-CPVMCL*(TCONV(1)-273.15) SX(1)=CPN(I)*TCONV(1)+GZ(1) HX(1)=SX(1)+LV(1)*RCONV(1) DO K = 2, NLEV CPN(K)=CPD*(1.-RCONV(K))+RCONV(K)*CPV TVX=TCONV(K)*(1.+RCONV(K)*EPSI-RCONV(K)) TVY=TCONV(K-1)*(1.+RCONV(K-1)*EPSI-RCONV(K-1)) GZ(K)=GZ(K-1)+0.5*RD*(TVX+TVY) : *(PPLAY(K-1)-PPLAY(K))/PAPRS(K) LV(K)=LV0-CPVMCL*(TCONV(K)-273.15) SX(K)=TCONV(K)*CPN(K)+GZ(K) HX(K)=SX(K)+LV(K)*RCONV(K) ENDDO DO K = 1, NLEV SXG1(J,K)=SX(K)/1000. HXG1(J,K)=HX(K)/1000. ENDDO CC (acz) interactive SST and STL: computed at each CC gridpoint (then masked to determine TS) CC with different QFLUX (defined in do_inphys) CC IF (.NOT. PRESC_SST) THEN SST_TEND = (-SH - LH + SSW(J) - SLW(J) + QFLUX(J)) : / HEATCAPOCEA SST1(J) = SST1(J) + SST_TEND * DELT END IF IF (.NOT. PRESC_STL) THEN STL_TEND = (-SH - LH + SSW(J) - SLW(J) + QFLUX(J)) : / HEATCAPLAND STL1(J) = STL1(J) + STL_TEND * DELT END IF CC CC (nprive+acz) Interactive soil moisture CC c convert latent heat flux to evaporation and c evaporation in kg m-2 s-1 into evaporation in cm/s EBKT = LH / (10. * ALV) c convert precipitation in mm/day in cm/s PBKT = (PRECIP + PRADJ)/(10*24*3600) IF (.NOT. PRESC_BKT) THEN BUCKET(J) = BUCKET(J) + (PBKT-EBKT) * DELT IF (BUCKET(J) .LT. 0) then BUCKET(J) = 0 ELSEIF (BUCKET(J) .GT. BKT0) then BUCKET(J) = BKT0 ENDIF ELSE BUCKET(J) = BKT0 ENDIF C ================================================================ 9999 CONTINUE ! loop over atmospheric columns ends up here C ================================================================ C DUMMY = 0 C DO k = 1, K_PBL C DUMMY = DUMMY + UT_PBL_W(100,k) + UT_PBL_E(100,k) C END DO C write(*,*) 'MITPHYS_DRIVER:' C write(*,*) 'UT_PBL(100,1:K_PBL):', DUMMY C write(*,*) 'USTR(100) = ',USTR(100) C DUMMY = 0 C DO k = 1, K_PBL C DUMMY = DUMMY + VT_PBL_S(100,k) + VT_PBL_N(100,k) C END DO C write(*,*) 'VT_PBL(100,1):', DUMMY C write(*,*) 'VSTR(100) = ', VSTR(100) C write(*,*) ROWS C C write(*,*) 'PHY_DRIVER end, TG1(65,1):', TG1(65,1) C write(*,*) 'PHY_DRIVER end, CBMFG1(85):', CBMFG1(85) FirstCall = .FALSE. RETURN END C $Id: mitphy_driver.F,v 1.1 2003/08/20 15:24:59 czaja Exp $