C ---------------------------------------------------------------------- C COUPLED ATMOSPHERIC BOUNDARY LAYER-PLANT-SOIL (CAPS) MODEL C ---------------------------------------------------------- C ONE-DIMENSIONAL ATMOSPHERIC BOUNDARY LAYER MODEL THAT SIMULATES C THE ATMOSPHERE, SOIL, AND VEGETATED SURFACE, VERSION 1.0.4C. C ---------------------------------------------------------------------- C REFERENCE: A USER'S GUIDE TO OSU1DPBL, VERSION 1.0.4, MAR 1991. C REFER TO THIS GUIDE FOR A COMPLETE DESCRIPTION OF SUBROUTINE C ARGUMENTS, COMMON BLOCK ELEMENTS, AND OTHER INFORMATION. C C ORIGINAL REFERENCE: TROEN AND MAHRT, 1986: "A SIMPLE MODEL OF C THE ATMOSPHERIC BOUNDARY LAYER; SENSITIVITY TO SURFACE C EVAPORATION." BOUNDARY LAYER METEOROLOGY, 37, 129-148. C ---------------------------------------------------------------------- C MODIFICATIONS FROM 1.0.1 (JUNE 1988): C ---------------------------------------------------------------------- C 02-AUG-88, P. RUSCHER - VERSION 1.0.1B C * DEFAULT GRID NOW NUMBER 1 (20 METER RESOLUTION). C * TURNED OFF CLOUD-RADIATION FEEDBACK IN SUBROUTINE SUN. C * CAP PLACED ON THE COUNTERGRADIENT CORRECTION FACTOR FOR HEAT C (CGH); COUNTERGRADIENT CORRECTION FOR MOISTURE ELIMINATED C (CGW = 0). C * CRITICAL RICHARDSON NUMBER (RICR) = 0.5. c Note (3/14/90 J. Kim) c 1. Grid #1 is extended to 9 km c 2. Default critical Richardson number of set back to 1, and can c be modified via interactive I/O prior to the main calculation. C ---------------------------------------------------------------------- C 29-NOV-88, WAYNE GIBSON, V - 1.0.2 C * MODIFIED TO RUN ON THE MACINTOSH. c Note (3/14/90 J. Kim) c This version is strictly for VAX or comparable machine. c all of the MacIntosh stuffs are removed from the code. C ---------------------------------------------------------------------- C 16-JAN-89, MIKE EK - V. 1.0.2 C * ALBEDO, TWILT, TSOFEF, TSO0, SIGMAF, PC, ISOIL, CMC, SCANOP, C CLC, IFCLD, AND IFTC INITIALIZED IN (INPUT) CONTROL FILE C RATHER THAN IN BLOCK DATA STATEMENTS. SEE USER'S GUIDE FOR C A COMPLETE DESCRIPTION OF THE CONTROL FILE. C ---------------------------------------------------------------------- C WINTER-89/90, M. EK - V. 1.0.2A TO 1.0.2G C * MISCELLANEOUS TESTS AND CHANGES TO COMMENTS. C ---------------------------------------------------------------------- C 26-JAN-89, M. EK - V. 1.0.2H C * CHANGED SOIL HEAT CAPACITY CODE IN SUBROUTINE HRT TO INCLUDE C HEAT CAPACITY OF AIR IN SOIL. C ---------------------------------------------------------------------- C WINTER-89/90, M. EK - V. 1.0.2I TO 1.0.2J C * MISCELLANEOUS TESTS AND CHANGES TO COMMENTS. C ---------------------------------------------------------------------- C 20-APR-89, M. EK - V. 1.0.2K C * REMOVED CALCULATION OF FRACTIONAL CLOUD COVER FROM C SUBROUTINES PBLL AND SUN. FRACTIONAL CLOUD COVER NOW C CALCULATED IN SUBROUTINES STBCLC (CLOUD COVER UNDER STABLE C CONDITIONS) AND UNSCLC (CLOUD COVER UNDER UNSTABLE C CONDITIONS), AND CALLED FROM PBLL. FRACTIONAL CLOUD COVER C (CLC) PLACED IN COMMON BLOCK 'RAD' FOR LATER USE IN C SUBROUTINE SUN. C ---------------------------------------------------------------------- C 15-MAY-89, M. EK - V. 1.0.2K C * CLOUD-RADIATION INTERACTION FLAG ADDED TO COMMON BLOCK RAD C FOR USE IN SUBROUTINE SUN. THIS ELIMINATES THE NEED TO SET C CLOUD COVER TO ZERO IN SUBROUTINE SUN IN ORDER TO TURN OFF C CLOUD- RADIATION INTERACTION (EARLIER CALLED CLOUD-RADIATION C FEEDBACK). C * MOST ALL COMMON BLOCK ELEMENTS ARE NOW IDENTICAL THROUGHOUT C OSU1DPBL. SEE USER'S GUIDE/MODEL CHANGES FOR REMAINING C LIST. C ---------------------------------------------------------------------- C 20-MAY-89, M. EK - V. 1.0.2L C * VERTICAL ADVECTION OF MOMENTUM (W*DU/DZ, W*DV/DZ) ADDED TO C THE MODEL IN THE PROGRAM DRIVER. C ---------------------------------------------------------------------- C SPRING-89, M. EK - V. 1.0.2M TO 1.0.2U C * MISCELLANEOUS TESTING AND CHANGES TO COMMENTS. C ---------------------------------------------------------------------- C 28-JUN-89, M. EK - V. 1.0.2V C * PARAMETERIZED RADIATION FLAG (IFPR) ELIMINATED FROM CONTROL C FILE AND REPLACED BY THE CLOUD-RADIATION INTERACTION FLAG C (IFCRI). THIS ALLOWS CLOUDS TO BE TURNED ON/OFF QUICKER. C * PARAMETERIZED RADIATION IS NOW HARD-WIRED AS ZERO (OFF) IN C THE MODEL, AND HAS BEEN REMOVED AS AN ARGUMENT IN SUBROUTINE C SUN AND ADDED TO THE COMMON BLOCK 'RAD'. C ---------------------------------------------------------------------- C 15-AUG-89, M. EK - V. 1.0.2W C * IN SUBROUTINE SSTEP PUT CAP ON SOIL MOISTURE SO THAT IT C CANNOT BECOME SUPERSATURATED. EXCESS MOISTURE IS RUNOFF. C * APPARENT SURFACE SOIL MOISTURE (THSFC) NOW CALCULATED IN C DIRECT SOIL EVAPORATION SUBROUTINE (SEDIR). THSFC IS C COMPARED TO THE AIR-DRY VALUE (TSO0) TO DETERMINE IF DIRECT C SOIL EVAPORATION IS POTENTIAL OR CONTROLLED BY UPWARD C MOISTURE FLUX IN THE SOIL. PREVIOUSLY THSFC WAS NOT C CALCULATED TO DETERMINE DIRECT SOIL EVAPORATION; IT WAS DONE C IMPLICITLY WITHOUT THSFC. C ---------------------------------------------------------------------- C 19-AUG-89, BERT HOLTSLAG - V. 1.0.2X C * ROUGHNESS LENGTHS FOR MOMENTUM AND HEAT ARE DIFFERENT (SEE C SFLX). C ---------------------------------------------------------------------- C 23-AUG-89, B. HOLTSLAG AND M. EK - V. 1.0.2X C * MOVED THE CALCULATION OF THE SURFACE EXCHANGE COEFFICIENTS C AND BULK RICHARDSON NUMBER TO SUBROUTINE SFCXCH. C ---------------------------------------------------------------------- C 24-AUG-89, B. HOLTSLAG AND M. EK - V. 1.0.2X C * DISPLACEMENT HEIGHT (ZD0) ADDED TO MODEL. C ---------------------------------------------------------------------- C SUMMER-89, M. EK AND B. HOLTSLAG - V. 1.0.3 TO 1.0.3D C * MANY TEST VERSIONS AND NUMEROUS COSMETIC CHANGES TO COMMENTS C AND PROGRAM STRUCTURE. C ---------------------------------------------------------------------- C 29-AUG-89, B. HOLTSLAG AND M. EK - V. 1.0.3E C FOLLOWING HOLTSLAG ET AL. (1989) THE FOLLOWING CHANGES WERE C MADE: C * SIMILARITY PROFILE FUNCTIONS MODIFIED FOR UNSTABLE CASE, C * ADDED 2 METER CALCULATION OF TEMPERATURE AND WIND MODIFIED C BY STABILITY, C * SURFACE EXCHANGE COEFFICIENTS CHANGED USING LOUIS ET AL. C (1982) FOR UNSTABLE CASE, AND PRANDTL NUMBER = 1 FOR C NEUTRAL CONDITIONS, C * SURFACE VIRTUAL TEMPERATURE (IN SFCXCH) IS CALCULATED WITH C SURFACE MOISTURE INSTEAD OF FIRST ATMOSPHERIC MODEL LAYER C VALUE, C * BULK RICHARDSON NUMBER USES REFERENCE TEMPERATURE AT FIRST C ATMOSPHERIC MODEL LEVEL INSTEAD OF SKIN TEMPERATURE. C ---------------------------------------------------------------------- C 29-AUG-89, B. HOLTSLAG AND M. EK - V. 1.0.3E C * IFPR ELIMINATED FROM SUBROUTINE SUN. SATTERLUND DOWNWARD C LONGWAVE RADIATION CALCULATION NOW HARDWIRED (IFPR=1). C ---------------------------------------------------------------------- C 30-AUG-89, M. EK - V. 1.0.3E C * Z0H, ZD0, TSOREF, AND PC NOW INCLUDED IN THE OUTPUT FILE. C ALL COMMON BLOCK VARIABLES READ INTO THE CONTROL FILE ARE C INITIALIZED WITH DUMMY VALUES IN THE BLOCK DATA A1-A4. THIS C ELIMINATES A PROBLEM WITH THE MACTRAN+ FORTRAN COMPILER. C ---------------------------------------------------------------------- C 06-SEP-89, B. HOLTSLAG - V. 1.0.3E C * IF ON INPUT PC<0 THEN CANOPY RESISTANCE (RC) IS -PC AND PC C IS THEN CALCULATED FROM RC VIA PENMAN-MONTEITH IN SFLX. C OTHERWISE PC IS JUST TAKEN AS SPECIFIED IN THE (INPUT) C CONTROL FILE. AT THIS TIME THIS CHANGE IS NOT USED FOR THE C CASE OF SNOW. C ---------------------------------------------------------------------- C 07-SEP-89, M. EK - V. 1.0.3E C * ADD PC TO HOURLY OUTPUT. IF RC IS CONSTANT, PC IS C CALCULATED IN SFLX. C * ADD RC TO COMMON BLOCK SOIL1 AND TO HOURLY OUTPUT. IF PC IS C CONSTANT, RC IS CALCULATED IN SFLX. C ---------------------------------------------------------------------- C 12-SEP-89, M. EK - V. 1.0.3E C * SOLAR (INCOMING SOLAR RADIATION BEFORE IT REACHES THE C SURFACE) ADDED TO COMMON BLOCK 'RAD' AND TO HOURLY OUTPUT. C ---------------------------------------------------------------------- C 19-SEP-89, M. EK - V. 1.0.3E C * CHANGED CONTROL (INPUT) FILE TO INCLUDE SEPARATE VALUES FOR C SOIL MOISTURE AT THE DIFFERENT LAYERS INSTEAD OF HAVING A C CONSTANT INITIAL PROFILE. CURRENTLY THIS IS ONLY VALID FOR C A MODEL WITH TWO SOIL LAYERS. C ---------------------------------------------------------------------- C 02-OCT-89, M. EK - V. 1.0.3E C * CHANGED SOLAR TO SOLARN IN COMMON BLOCK 'RAD' AND IN HOURLY C OUTPUT, ALTHOUGH IT IS STILL CALLED 'SOLAR' IN HOURLY C OUTPUT. SOLARN IS THE NET SOLAR RADIATION AT THE SURFACE, C SOLARN = SOLAR * (1-ALBEDO). C ---------------------------------------------------------------------- C 05-OCT-89, M. EK - V. 1.0.3E C * BOTH CLOUD COVER MODELS ARE NOW USED REGARDLESS OF C STABILITY. THE FRACTIONAL CLOUD COVER IS JUST THE MAXIMUM C OF THE TWO. C * IFCRI FLAG CHANGED TO INCLUDE A VALUE OF -1. WHEN IFCRI=-1 C CLOUD COVER IS DIAGNOSED BUT THERE IS NO MODIFICATION OF C RADIATION IN SUBROUTINE SUN. (FOR IFCRI=0 IT IS THE SAME AS C IFCRI=-1 EXCEPT THAT CLOUD COVER IS NOT DIAGNOSED; FOR C IFCRI=1 CLOUD COVER IS DIAGNOSED AND THERE IS MODIFICATION C OF RADIATION.) C * WHEN A DISPLACEMENT HEIGHT (ZDO) IS SPECIFIED, DATA FROM THE C CONTROL FILE IS NOT USED FOR Z>0 AND Z<=ZDO. C ---------------------------------------------------------------------- C 14-FEB-90, M. EK - V. 1.0.3E C * MODIFICATIONS TO SUBROUTINE SFCXCH TO EXACTLY AND MORE C EXPLICITLY FOLLOW FORMULATIONS FROM LOUIS ET AL (1982) FOR C SURFACE EXCHANGE COEFFICIENTS. C * NOW ONLY USING SUBROUTINE STBCLC ROUTINE TO DETERMINE C FRACTIONAL CLOUD COVER. c --------------------------------------------------------------------- * 1 March 1990 by Jinwon Kim: Version := OSU2.FOR c * * Subroutine ROSR12 (tri-diagonal matrix inverter) now use * double-precision variables for internal calculation. * This change produces virtually same result as the previous * version (OSU1DPBL1_0_3F.FOR). * * This version is the starting point for the free-atmospheric * diffusion and gravity wave model c --------------------------------------------------------------------- * 2 March 1990 by Jinwon Kim: Version := OSU2_0.FOR c * * Structural changes in the subroutine "PBLL" * 1. Simplify * 2. Use "ROSR12" to invert tri-diagonal matrix. * Change of matrix solver can alter the result, but not * significantly (new method is probably better, and * necessary for the implication of free atmospheric * diffusion calculation. Previous matrix solver has * spurious noise in the free atmosphere, but with new * method, noise are suppressed significantly). c --------------------------------------------------------------------- * 12 March 1990 by J. Kim: Version := OSU2_1.FOR * * Critical Richardson number can be chosen as interactive I/O. * (default Ricr = 1.0) * * Flag for free atmospheric diffusion (IFREE) is installed, * and can be set interactively. * * When free atmospheric diffusion is not calculated, any * nonzero flux convergence above the PBL top is set * to be zero (they are noise of the matrix inversion). * * #1 grid is now extended to 9 km level. c Note --- * In the presence of prescribed mean vertical motion, which simulates * vertical advection by subsidence, upper levels near the domain top * becomes convectively unstable for a shallow domain (Ztop = 3 km). * Subsidence (imposed vertical motion) is necessary to prevent too * deep daytime PBL. c ---------------------------------------------------------------------- * 20 March 1990 Jinwon Kim Version := OSU3.FOR * * Parameterization for gravity wave momentum transfer is added * * Flag to turn on the wave drag is IGW=1. Otherwise, gravity * wave drag will be skipped. IGW also determines the type of * gravity wave model afterwards (by interactive I/O) * * For wave drag parameterization, see subroutine gw. * * Time scale to achieve the equailibrium state between the * wave breaking and mean flow adjustment (GWTIM) is built in. * By changing GWTIM one can reduce the forcing by gravity * wave drag (if GWTIM is less than time step). For detailed * procedure to prescribe the time scale, closely follow the * interactive I/O procedure. * * Common block added: * /gwave/vartpo,hktpo,gwtim,dudtw(nlevd),dvdtw(nlevd),tau(nlevd) c ---------------------------------------------------------------------- * 23 March 1990. Jinwon Kim. * * Four different versions of gravity wave models are in menu. * 1. Saturation/Linear LBC (subroutine gwls) * 2. Saturation/First-order LBC (subroutine gws) * 3. Supersaturation/Linear LBC (subroutine gwlss) * 4. Supersaturation/First-order LBC (subroutine gwss) c ---------------------------------------------------------------------- * 26 March 1990 J. Kim * * Height dependent/time independent value of geostrophic * wind can be assigned by setting input UGI(1) and VGI(1) * to be zero (both must be zero). Then, input u and v wind * profiles will be assigned as the vertical profile of * geostrophic wind. This is to minimize inertial oscillation. c ---------------------------------------------------------------------- * 25 July 1990 J. Kim * * Large-scale vertical motion can be given in two different * units; one is the old one (mb/hour) and the other is "m/s". * To use "m/s" unit for input "w", set "iwunit"=1 in the * control file. Otherwise, old unit will be assumed. c ..................................................................... c Nov 5 1990: Jinwon Kim for v. "osugw_105.sub" c * extra-fine grid (10 m spacing below 100 m level) in the low c level for finer resolution of very stable boundary layer. C ---------------------------------------------------------------------- C SPRING 90, M. EK - V. 1.0.3F - 1.0.3H C * TESTING OF NEW BOUNDARY-LAYER CLOUD COVER SCHEME. C ---------------------------------------------------------------------- C 13-AUG-90, M. EK - V. 1.0.3I C * TURNED OFF SUBROUTINE SSHEAT SO RH MAY EXCEED 1.O SO THAT C SPECIFIC HUMIDITY IS A TOTAL WATER SPECIFIC HUMIDITY; C NECESSARY FOR USE IN THE CLOUD COVER SCHEME. C ---------------------------------------------------------------------- C AUTUMN 90, M. EK - V. 1.0.3J - 1.0.3M C * TESTING OF NEW BOUNDARY-LAYER CLOUD COVER SCHEME. C ---------------------------------------------------------------------- C 14-FEB-91, M. EK - V. 1.0.4 C * CORRECTIONS IN SUBROUTINE SUN FOR DOWNWARD LONGWAVE C RADIATION IN THE PRESENCE OF CLOUDS. C * TRANSMISSION FUNCTION FOR SOLAR RADIATION THROUGH CLOUDS IN C SUBROUTINE FUNCTION FOLLOWS LIOU (1976) FOR SHALLOW CUMULUS, C DESCRIBED IN EK AND MAHRT (1991). C * NEW BOUNDARY-LAYER CLOUD COVER SCHEME VALID FOR BOTH STABLE C AND UNSTABLE REGIMES FOLLOWING EK AND MAHRT (1991). C ---------------------------------------------------------------------- C 20-MAR-91, J. KIM - V. 1.0.4A C * MODIFICATIONS TO THE SUBROUTINE PRTFLX USED FOR OUTPUT FROM C GRAVITY WAVE MODEL TESTING. C 20-Mar-91, M. EK C * MINOR CHANGES IN SUBROUTINE CLOUD, CHANGE 'MAX' TO 'IMAX' C (COMPILER SPECIFIC FOR THE NEXT WORKSTATION). C * CHANGE I/O SO THAT IT CAN RUN ON THE NEXT USING ABSOFT FTN. C ---------------------------------------------------------------------- C 27-JAN-92, M. EK - V. 1.0.4B C * REMOVE SEVERAL FILES FROM BEING WRITTEN (USED IN TESTING THE C GRAVITY WAVE PROGRAM). TURN OFF THE GRAVITY WAVE PROGRAM. C * CHANGE I/O SO THAT IT CAN RUN ON THE NEXT USING ABSOFT FTN. C ---------------------------------------------------------------------- C 28-MAY-92, M. EK - V. 1.0.4B1 C * REMOVE GRAVITY WAVE-RELATED SUBROUTINES IN THE MODEL. C 26-OCT-92, M. EK C * CHANGE CLOUD COVER FLAG. IFCRI=-1, CLOUD COVER FIXED C THROUGHOUT SIMULATION AT INITIAL VALUE IN CONTROL FILE. C ---------------------------------------------------------------------- C 10-DEC-92, M. EK and B. HOLTSLAG - V. 1.0.4B1_KNMI - side test C * VARIABLE CANOPY RESISTANCE USING SUBROUTINE CANRES.F CALLED C FROM SUBROUTINE SFLX. C * ADD FILE 'RCTAB' AS OUTPUT WITH CANOPY RESISTANCE TERMS. C * RICR = 0.25; appropriate when using fine vertical grid. C ---------------------------------------------------------------------- C 12-MAY-93, M. EK and B. HOLTSLAG - V. 1.0.4B1_KNMI2 - side test C * TURN ON THE COUNTERGRADIENT TERM FOR MOISTURE (CGW) C ---------------------------------------------------------------------- C 13-MAY-93, M. EK and B. HOLTSLAG - V. 1.0.4B1_KNMI2a - side test C * use 2nd model level virtual potential temperature instead c of virtual temperature in boundary-layer height calc. c * for reference pressure, define it equal to the surface c pressure rather than 1000 hPa. Important!!! - later, c use reference pressure=1000 hPa which requires changing c surface energy balance/Tsfc, Theta-sfc portions of model c code. C ---------------------------------------------------------------------- C 18-MAY-93, M. EK - V. 1.0.4B2_KNMI C * FIX POTENTIAL EVAPORATION AND SURFACE TEMPERATURE FORMULAS C TO INCLUDE A 1000 MB-BASED POTENTIAL TEMPERATURE. C * MISC. FIXES IN SOIL WATER RELATIONS AND SUBSURFACE RUNOFF; C COMPARE WITH VERSION 1.0.4B1. C * FIX CANOPY RESISTANCE/PLANT COEFFICIENT RELATIONSHIP DUE C TO THE FIX TO SURFACE POTENTIAL TEMPERATURE. C ---------------------------------------------------------------------- C 05-NOV-93, M. EK - V. 1.0.4B2_KNMI C * FIX MORE CODE RELATED TO T-SFC NOT EQUAL TO THETA-SFC C ---------------------------------------------------------------------- C 08-NOV-95, M. EK - V. CAPS_104C C * FIX SOME CODE RELATED TO CALCULATION OF TH2 (POTENTIAL TEMP C AT FIRST ATMOS. MODEL LEVEL; DEFAULT=20M). C ('CAPS' model = Coupled Atmospheric Boundary Layer-Plant-Soil model) C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT OR INITIALIZED: C /FIELDS/ UNM, VNM, RNM, QNM, UG, VG C /GRID/SS, SSK (CURRENT), NOOFGR, MZ, MZM C /SOIL/ WSOIL, CMC, NSOIL, TSOIL, ESD C /SOIL1/ IFTC, TSO0, SIGMAF, TSOSAT, TWILT, SCANOP, PC, TSOREF C /SOIL2/ DSOIL, ISOIL C /PBL/ ZZ (CURRENT), WINIT, RICR, PINK, DTNW, KOOL, IFCLD C /SFCL/ Z0, Z0H, ZDO, TSNOW C /AUXIL/ PSFC, PRCP, UFLX, VFLX, TFLX, QFLX C /RAD/ FDOWN, MO, DY, CLC, TREF, SLO, SLA, ALBEDO, IFCRI C /AUX12/ ETOT C /XLUN/ LUNC C /HCALC/ TIMEIS (CURRENT) c ----------------------------------------------------------------------- C OUTPUT OR MODIFIED: C /FIELDS/ UN, UNM, VN, VNM, RN, RNM, QN, QNM C /GRID/ SSK (UPDATED) C /SOIL/ SSOIL C /PBL/ DUDT, DVDT, DTDT, DWDT, ZZ (UPDATED), MZBL, DELT2 C /SFCL/ TS, WS, U2, V2, T2, TH2, W2, ZZ2 C /AUXIL/ DEW, ACCP, ACCD C /RAD/ TSUN C /XLUN/ LUNB, LUND C /HCALC/ TIMEIS (UPDATED) C ---------------------------------------------------------------------- C BLOCK DATA USAGE: C A1 RADIATION PARAMETERS C A2 SOIL MODEL PARAMETERS C A3 COMPUTATION TABLE FOR SATURATION VAPOR PRESSURE C A4 BOUNDARY LAYER MODEL AND I/O PARAMETERS C ---------------------------------------------------------------------- C LOOP STRUCTURES FOR MAIN PROGRAM DRIVER (PBLDRV): C 1xx Input C 2xx Initialization C 3xx Time Stepping C 4xx Tendency Computations C ... C 8xx Output formats C 9xx Flow control C ---------------------------------------------------------------------- C SUBROUTINES CALLED: WINTB (FOR MACINTOSH), TIME AND DATE (FOR C VAX), ZTOSIG, HFAK, CHACHA, SIGTOZ, SUN, C SFLX, PBLL, SVP, SSHEAT, ERR1, PRINT C ********************************************************************** PROGRAM PBLDRV C ---------------------------------------------------------------------- C NLEVD=LEVELS, NLEVI=INPUT LEVELS, NSOLD=SOIL LEVELS. C ---------------------------------------------------------------------- parameter(NLEVD=100,NLEVI=100,NSOLD=10,ngrid=8) COMMON/FIELDS/UN(NLEVD),UNM(NLEVI),VN(NLEVD),VNM(NLEVI), . RN(NLEVD),RNM(NLEVI),QN(NLEVD),QNM(NLEVI), . UG(NLEVD),VG(NLEVD) COMMON/GRID/SS(NLEVD),SSK(NLEVD),MZ,MZM,NOOFGR COMMON/SOIL/WSOIL(NSOLD),CMC,NSOIL,SSOIL,TSOIL(NSOLD),ESD COMMON/SOIL1/IFTC,TSO0,SIGMAF,TSOSAT,TWILT,SCANOP, . PC,RC,KWILT,CFACTR,TSOREF COMMON/SOIL2/DSOIL(NSOLD),ISOIL,IFSOIL COMMON/PBL/DUDT(NLEVD),DVDT(NLEVD),DTDT(NLEVD),DWDT(NLEVD), . VERT(NLEVD),ZZ(NLEVD),WINIT(NLEVD),HPBL,MZBL,DELT2, . RICR,PINK,DTNW,KOOL,PBLK(NLEVD),CLOUDK(NLEVD), . ZLCL,CLTOP,DCZ,IFCLD COMMON/SFCL/TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,TSNOW,Z01,CBAGK COMMON/AUXIL/PSFC,HEAT,TAOX,TAOY,EVAP,XL,CG,THERM,DELTAT,FH, . PRCP,DEW,ACCP,ACCD,UFLX(NLEVD),VFLX(NLEVD), . TFLX(NLEVD),QFLX(NLEVD) COMMON/RAD/FDOWN,MO,DY,TSUN,CLC,TREF,SLO,SLA,ALBEDO,IFCRI,IFPR, . SOLARN COMMON/SOIL3/B(11),SATPSI(11),SATKT(11),TSAT(11) COMMON/AUXI2/EP,ETOT,PR COMMON/XLUN/LUNC,LUNB,LUND COMMON/HCALC/MZH,TIMEIS character FINTC*64,FOUT*64,FOUTD*64,FNAME*64,FEXT*64 character tab*1 character*72 texti dimension dz(nlevi),ugi(nlevd),vgi(nlevd),tgi(nlevd),tg(nlevd), . vertm(nlevi) dimension inli(ngrid),igrid(nlevd,ngrid),dqdz(nlevd),dtdz(nlevd), . dudz(nlevd),dvdz(nlevd) real envk(nlevd),work(nlevd),p(nlevd) logical kwilt c -------------------------------------------------- c Set up data for different grid in units of ZUNIT c INLI is the number of grid points for each grid. c c #1 grid is extended upto 10 km level c #7 grid is a coarse one extending upto 1 km. c 3/26/1990 J. Kim c -------------------------------------------------- data inli/70,25,13,7,4,9,27,75/ data igrid/ . 0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40, . 45,50,55,60,65,70,75,80,85,90,95,100,105,110,115,120,125, . 130,135,140,150,160,170,180,190,200,220,240,260,280,300,320, . 340,360,380,400,430,460,500,550,600,650,700,750,800,850,900, . 950,1000,30*0, . 0,1,2,3,4,5,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40, . 80,75*0, . 0,1,2,4,7,10,14,20,26,32,40,60,80,87*0, . 0,1,5,10,20,40,80,93*0, . 0,1,40,80,96*0, . 0,4,16,35,71,132,213,312,430,91*0, . 0,4,10,20,30,45,60,85,100,150,200,250,300, . 350,400,450,500,550,600,650,700,750,800,850,900, . 950,1000,73*0, . 0,1,2,3,4,5,6,7,8,9,10,12,14,16,18,20,22,24,26,28,30,32,34,36, . 38,40,45,50,55,60,65,70,75,80,85,90,95,100,105,110,115,120,125, . 130,135,140,150,160,170,180,190,200,220,240,260,280,300,320, . 340,360,380,400,430,460,500,550,600,650,700,750,800,850,900, . 950,1000,25*0/ data zunit/10./ data init,mz0,xskip,idgeo/0,1,777.,0/ data grd,mbtoka,g,cp/0.034146341,100,9.806,1004.5/ data rkfac/1./ c ........................................................... pi=acos(-1.) pi2=2.*pi tab = char(9) c ----------------------------------------------------------- c Height at grid level for #1 and #6 grid is IGRID*10 m. c For grids #2 - #5, heights at grid levels are IGRID*50 m. c ----------------------------------------------------------- do100 k=1,nlevi do100 kg=2,5 igrid(k,kg)=igrid(k,kg)*5 100 continue c -------------------------------- c Set parameters for evaporation c -------------------------------- DTNW=86400. ETOT=0. EVAP=0. ACCP=0. ACCD=0. c Free atmospheric diffusion flag = IFREE ifree = 1 ifkpbl = 0 c ---------------------------------------------------------------- write(*,*) write(*,*)' PBL MODEL PROGRAM RUN ' c ---------------------------------------------------- * This block will run under VAX/VMS operating system c ---------------------------------------------------- WRITE(*,*) WRITE(*,*) WRITE(*,800)' CONTROL FILE NAME --> ' READ(*,810) FINTC WRITE(*,800)' OUTPUT FILE NAME ---> ' READ(*,810) FOUT LENF = LEN(FINTC) KSKIP = 0 KFB = 0 DO110 K=1,LENF IF (FINTC(K:K).EQ.']') KFB = K IF (FINTC(K:K).EQ.' '.AND.K.NE.1.AND.KSKIP.EQ.0) THEN KFE = K-KFB-1 KSKIP = 1 ENDIF 110 CONTINUE LENF = KFE FNAME = FINTC(KFB+1:LENF-4) LENF = LENF-4 FEXT(1:4) = '.dmp' FOUTD(1:LENF) = FNAME(1:LENF) FOUTD(LENF+1:LENF+4) = FEXT(1:4) LENF = LENF+4 WRITE(*,*) WRITE(*,*)' An initial dump of all the controlling parameters' WRITE(*,*)' is made to file ',FNAME(1:LENF),' and can be ' WRITE(*,*)' immediately viewed to check input parameters.' WRITE(*,*) WRITE(*,*)'program executing' WRITE(*,*) c ------------------------------------ * open data files for input/output - c ------------------------------------ open(unit=lunc,file=fintc,status='old',form='formatted') cx open(unit=lunb,file=fout,status='new',form='unformatted') open(unit=lunb,file=fout,status='unknown',form='unformatted') cx open(unit=lunb,file=fout,status='unknown',form='unformatted', cx . block=-1) open(unit=lund,file=foutd,status='unknown',form='formatted') open(unit=12,file='rctab',status='unknown',form='formatted') ircflg = 0 iflag = 0 write (12,*) '*' write (12,*) 'time',tab,'rcs',tab,'rct',tab,'rcq',tab, . 'rcsoil',tab,'rc',tab,'PC' c -------------------------------------------------- * Read interactive input data and dump the final - * parameters in the input-data dump file. - c -------------------------------------------------- read(lunc,*,end=990) aaa READ(LUNC,*,END=990) IFHF,IFCRI,IFSNO,IFCLD WRITE(LUNB) IFHF,IFCRI,IFSNO,IFCLD WRITE(LUND,*) 'PROGRAM OSUPBL1D DUMP OF INITIAL DATA' WRITE(LUND,*) WRITE(LUND,900) IFHF,IFCRI,IFSNO,IFCLD IF((IFHF .EQ. 0) .AND. (IFSNO .EQ. 1)) THEN WRITE(*,*)' INVALID COMBINATION OF FLAGS !!!!' WRITE(*,*)' NO HEAT FLUX WITH SNOW.' WRITE(*,*)' PROGRAM TERMINATED !!!!!!!!' STOP ENDIF READ(LUNC,*,END=990) DELTAT,TEND MZ=INLI(NOOFGR) WRITE(LUNB) NOOFGR,DELTAT,TEND,RICR,PINK,KOOL WRITE(LUND,905) NOOFGR,DELTAT,TEND,RICR,PINK,KOOL c ------------------------------------------------------------------- * Wind and soil modifications factors not supported in this version. c ------------------------------------------------------------------- READ(LUNC,840)TEXTI READ(LUNC,*,END=990)SLA,SLO,TZONE,Z0,Z0H,ZD0,ALBEDO WRITE(LUNB)TEXTI WRITE(LUND,910)TEXTI c ---------------------------------------------------------- * Modify default albedo if 1d model is run over snow/water. c ---------------------------------------------------------- IF(IFSNO.EQ.1)THEN ALBEDO=0.7 ENDIF IF(Z0.EQ.0.)THEN ALBEDO=0.1 ENDIF WRITE(LUNB)SLA,SLO,TZONE,Z0,Z0H,ZD0,ALBEDO WRITE(LUND,915)SLA,SLO,TZONE,Z0,Z0H,ZD0,ALBEDO READ(LUNC,*,END=990)MO,DY,TIMEIS WRITE(LUNB)MO,DY,TIMEIS WRITE(LUND,920)MO,DY,TIMEIS READ(LUNC,*,END=990)PSFC WRITE(LUNB)PSFC WRITE(LUNB)TREF READ(LUNC,*,END=990)CLC WRITE(LUNB)CLC READ(LUNC,*,END=990)CMC,SCANOP WRITE(LUNB)CMC write(lund,925)psfc,tref,clc,cmc READ(LUNC,*,END=990)PRST,PREND,PRCIP,ESD WRITE(LUNB)PRST,PREND,PRCIP,ESD,TSNOW WRITE(LUND,930)PRST,PREND,PRCIP,ESD,TSNOW READ(LUNC,*,END=990)ISOIL,TWILT,SIGMAF,IFTC,TSO0,TSOREF,PC WRITE(LUNB)ISOIL,TWILT,SIGMAF,IFTC,TSO0,TSOREF,PC,RC WRITE(LUND,935)ISOIL,TWILT,SIGMAF,IFTC,TSO0,TSOREF,PC READ(LUNC,*,END=990)WSOIL(1),WSOIL(2) READ(LUNC,*,END=990)NUG WRITE(LUND,940)NUG DO 130 I=1,NUG READ(LUNC,*,END=990)UGI(I),VGI(I),TGI(I) WRITE(LUND,945)UGI(I),VGI(I),TGI(I) 130 CONTINUE READ(LUNC,*,END=990)mz1,iwunit c ------------------------------------------------------------ * Read-in the profile data: z,u,v,w,t,q: Change to MKS units. * Allow for the displacement height as in zz(j) below. * If there is data in the control file between the surface * (z>0) and the displacement height (z<=zd0), then throw it * out by reading the next level of data in its place. c ------------------------------------------------------------ DO150 j=1,mz1 145 READ(LUNC,*,END=990)DZ(J),UNM(J),VNM(J),VERTM(J),RNM(J), . QNM(J) WRITE(LUND,955)DZ(J),UNM(J),VNM(J),VERTM(J),RNM(J),QNM(J) RNM(J)=RNM(J)+273.16 QNM(J)=QNM(J)*0.001 DZ(J)=DZ(J)-ZD0 IF(DZ(J).LT.0.)THEN DZ(J)=0. IF(J.GT.1)THEN mz1=mz1-1 GOTO 145 ENDIF ENDIF 150 CONTINUE WRITE(LUND,950)mz1 c ------------------------------------------------------- * Define soil temperatureas equal to bottom temperature, * and initialize the distribution of soil moisture. c ------------------------------------------------------- WRITE (LUND,960) NSOIL DO 160 JS=1,NSOIL TSOIL(JS) = RNM(1) WRITE (LUND,965) DSOIL(JS),WSOIL(JS),TSOIL(JS) 160 CONTINUE WRITE (LUNB) NUG DO 170 I=1,NUG WRITE (LUNB) UGI(I),VGI(I),TGI(I) 170 CONTINUE c ------------------------- * Read-in soil parameters c ------------------------- WRITE (LUNB) NSOIL DO 180 JS=1,NSOIL WRITE (LUNB) DSOIL(JS),WSOIL(JS),TSOIL(JS) 180 CONTINUE WRITE (LUNB) mz1 texti(71:72) = '1V' C ---------------------------------------------------------------------- C CLOSE FIXED DATA FILE,CLOSE FORMATTED FILE FOR DUMPING THE DATA, C OPEN THE BINARY OUTPUT FILE FOR MODEL RESULTS. C ---------------------------------------------------------------------- WRITE (LUND,820) FINTC WRITE (LUND,830) FOUT c ................................. c close input and dump data files c ................................. close(unit=lunc) close(unit=lund,status='keep') C ---------------------------------------------------------------------- C COMPUTE GRID LEVELS AND SUBTRACT DISPLACEMENT HEIGHT FOR C VEGETATION. C ---------------------------------------------------------------------- do200 j=1,mz zz(j)=zunit*float(igrid(j,noofgr))-zd0 if(zz(j).lt.0.)zz(j)=0. 200 continue MM=MZ DO 210 J=2,MZ IF ((DZ(MZ1) .LE. ZZ(J)) .AND. (DZ(MZ1) .GE. ZZ(J-1))) THEN MM=J GOTO 215 ENDIF 210 CONTINUE 215 IF (DZ(MZ1) .LT. ZZ(MM)) THEN MZ = MM - 1 ELSEIF (DZ(MZ1) .EQ. ZZ(MM)) THEN MZ = MM ENDIF WRITE (LUNB) MZ WRITE (*,*) ' MODEL LEVEL =', MZ c --------------------------- c CORIOLIS PARAMETER. - c --------------------------- FH=(1.458E-4)*SIN(0.0175*SLA) c --------------------------------------------------------------------- c END TIME IS CURRENT TIME PLUS SPECIFIED INTEGRATION PERIOD. - c --------------------------------------------------------------------- TEND = TEND + TIMEIS C ---------------------- C CHANGE GRID. - C ---------------------- DO220 J=1,MZ UN(J) = XSKIP VN(J) = XSKIP RN(J) = XSKIP QN(J) = XSKIP VERT(J) = XSKIP 220 CONTINUE UN(1) = 0. VN(1) = 0. C ---------------------------------------------------------------------- C IF IFIT=1 THEN NEW GRID FITTED BY LEAST SQUARES; C IF IFIT=0 THEN NEW GRID FITTED BY LINEAR INTERPOLATION. C ---------------------------------------------------------------------- IFIT=0 if(ifit.le.0)then do230 j=1,mz CALL CHACHA(UN,ZZ,J,J,UNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(VN,ZZ,J,J,VNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(RN,ZZ,J,J,RNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(QN,ZZ,J,J,QNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(VERT,ZZ,J,J,VERTM,DZ,MZ0,MZ1,XSKIP,UG) 230 continue else CALL CHACHA(UN,ZZ,MZ0,MZ,UNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(VN,ZZ,MZ0,MZ,VNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(RN,ZZ,MZ0,MZ,RNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(QN,ZZ,MZ0,MZ,QNM,DZ,MZ0,MZ1,XSKIP,UG) CALL CHACHA(VERT,ZZ,MZ0,MZ,VERTM,DZ,MZ0,MZ1,XSKIP,UG) endif C ------------------------------------ C EXTEND MZ AND DEFINE MZBL. - C ------------------------------------ NEWMZ = MZ IF(DZ(MZ1) .GT. ZZ(MZ)) THEN DO240 K = 1, MZ1 IF(DZ(K) .GT. ZZ(MZ)) THEN NEWMZ = NEWMZ + 1 IF(NEWMZ .GT. NLEVD) STOP ' TOO MANY INPUT LEVELS' ZZ(NEWMZ) = DZ(K) UN(NEWMZ) = UNM(K) VN(NEWMZ) = VNM(K) RN(NEWMZ) = RNM(K) QN(NEWMZ) = QNM(K) VERT(NEWMZ) = VERTM(K) ENDIF 240 CONTINUE ENDIF MZBL = MZ MZ = NEWMZ MZM = MZBL - 1 C ---------------------------------------------------------------------- C CHANGE TO SIGMA LEVELS, COMPUTE FACTORS IN HYDROSTATIC EQUATION. C ---------------------------------------------------------------------- CALL ZTOSIG(SS,RN,QN,ZZ,MZ) CALL HFAK(SS,MZ) C ---------------------------------------------------------------------- C ZERO SETTINGS AND COMPUTE (P/P*)**(-KAPPA) FOR THETA COMPUTATION. C ---------------------------------------------------------------------- DUDT(1) = 0. DVDT(1) = 0. DTDT(1) = 0. DWDT(1) = 0. DO250 J=1,MZ UFLX(J) = 0. VFLX(J) = 0. TFLX(J) = 0. QFLX(J) = 0. WINIT(J) = QN(J) p(j)=PSFC * SS(J) SSK(J) = (p(j)/1.E5)**(-0.286) vert(j) = vert(j)*aaa if(iwunit.ne.1)then vert(j)=-vert(j)*rn(j)*mbtoka/(pres*grd*3600.) endif 250 CONTINUE c --------------------------------------------------------------- c Current hour is stored to control the print outs (once/hour). c --------------------------------------------------------------- ih2=ifix(timeis) c --------------------------------------------------------------- c Setup for the first time step c ----------------------------- c Use small step for the first forward (Euler) scheme. c c Time step DELT2 enters through common /PBL/ into implicit c scheme in PBL model. Time step in soil model is DELTAT. c c # of levels in PBL is MZBL which may be different from total c model levels (MZBL = MZ or MZ-1). c 10 MAY 1988 P. RUSCHER c --------------------------------------------------------------- c All turbulent fluxes are computed upto MZH if free atmospheric c diffusion is turned off (IFREE.ne.1). If free atmospheric c diffusion is turned on (IFREE = 1), turbulent fluxes are c calculated upto the layer MZBL-1. c c Gravity wave drag is calculated within the layer between the c prescribed mountain top (vartpo) and the domain top or the c critical (or unstable) layer if critical (or unstable layer) c is present between the mountop and domain top levels. c 22 March 1990 Jinwon Kim c --------------------------------------------------------------- 301 deltat=deltat/4. delt2=deltat soildt=deltat c -------------------------------------------- c Compute current value of geostrophic wind. c -------------------------------------------- tg(1)=timeis ug(1)=xskip vg(1)=xskip ns1=1 ns2=1 call chacha(ug,tg,ns1,ns2,ugi,tgi,mz0,nug,xskip,dz) if(ns1.ne.1.and.ns2.ne.1)call err1(ns1,ns2) call chacha(vg,tg,ns1,ns2,vgi,tgi,mz0,nug,xskip,dz) if(ns1.ne.1.and.ns2.ne.1)call err1(ns1,ns2) c ---------------------------------------------- c If in precipitation period, turn on the rain, c and accumulate rainfall. c ---------------------------------------------- IF ((TIMEIS .GE. PRST) .AND. (TIMEIS .LT. PREND)) THEN PRCP=PRCIP else prcp=0. ENDIF ACCP = ACCP + PRCP*DELTAT C ------------------------ C Update model time - C ------------------------ timeis=timeis+deltat/3600. C ---------------------------------------------------------- C SET TWO TIMELEVELS EQUAL BEFORE MAKING FORWARD STEP - C from unm to un, etc. (Computational I.C.) - C ---------------------------------------------------------- DO 310 IJ=1,MZ UNM(IJ) = UN(IJ) VNM(IJ) = VN(IJ) RNM(IJ) = RN(IJ) QNM(IJ) = QN(IJ) VERTM(IJ) = VERT(IJ) 310 CONTINUE C ------------------------------------------------------------------ C SOLVE HYDROSTATIC EQUATION TO OBTAIN Z AT SIGMA LEVELS. - C SET SFC VALUES AND VALUES AT FIRST LEVEL FOR USE IN SFLX. - C ------------------------------------------------------------------ CALL SIGTOZ(ZZ,RN,QN,MZ) TS = RNM(1) T2 = RNM(2) C ---------------------------------------------------------- C TH2 IS THE DRY-ADIABATIC TEMPERATURE FOR A PARCEL - C descending from zz2 to surface. - C ---------------------------------------------------------- cx TH2 = T2 * SS(2)**(-0.286) TH2 = T2 * SSK(2) WS = QNM(1) W2 = QNM(2) U2 = UNM(2) V2 = VNM(2) ZZ2 = ZZ(2) C ---------------------------------------------------------------------- C COMPUTE DOWNWARD RADIATION. C COMPUTE SURFACE DRAG COEFFICIENTS, UPDATE SOIL AND SFC TEMP, C SFC MOIST. ACCUMULATE DEWFALL. C ---------------------------------------------------------------------- TSUN = TIMEIS + TZONE CALL SUN c CALL SFLX (FDOWN,SSOIL,PSFC,PRCP,DEW,SOILDT,IFHF,IFSNO,ircflg) CALL SFLX . (FDOWN,SSOIL,PSFC,PRCP,DEW,SOILDT,IFHF,IFSNO,ircflg,iflag) ACCD = ACCD + DEW*DELTAT C ---------------------------------------------------------------------- C LOAD PBL COLUMN INTO D( )DT WORKARRAYS, CALL PBL MODEL ON EXIT C D( )DT ARRAYS CONTAIN TENDENCIES DUE TO PBL MIXING. C ---------------------------------------------------------------------- QNM(1) = WS QN(1) = WS DO400 J=2,MZBL jp=j+1 DZZ = ZZ(jp) - ZZ(J) DUDT(J) = UNM(J) DVDT(J) = VNM(J) DTDT(J) = RNM(J)*SSK(J) DWDT(J) = QNM(J) IF (J .LT. MZBL) THEN DUDZ(J) = (UNM(jp)-UNM(J))/DZZ DVDZ(J) = (VNM(jp)-VNM(J))/DZZ DQDZ(J) = (QNM(jp)-QNM(J))/DZZ DTDZ(J) = (RNM(jp)-RNM(J))/DZZ + G/CP ENDIF 400 CONTINUE CALL PBLL(ifree,envk,ifkpbl,rkfac,rb2) c ------------------------- * DO FORWARD STEP. - c ------------------------- RN(1) = TS QN(1) = WS do320 ij=2,mzbl UGZ=UG(1) VGZ=VG(1) uadv=vert(ij)*dudz(ij) vadv=vert(ij)*dvdz(ij) tadv=vert(ij)*dtdz(ij) qadv=vert(ij)*dqdz(ij) fv=fh*(vnm(ij)-vgz) fu=fh*(ugz-unm(ij)) UN(IJ)=UNM(IJ)+(fv-uadv+DUDT(IJ))*DELTAT VN(IJ)=VNM(IJ)+(fu-vadv+DVDT(IJ))*DELTAT RN(IJ)=RNM(IJ)+(DTDT(IJ)/SSK(IJ)-tadv)*DELTAT QN(IJ)=QNM(IJ)+(DWDT(IJ)-qadv)*DELTAT c CALL SVP(QSIJ,ESIJ,p(ij),RN(IJ)) QN(IJ) = AMAX1(0.,QN(IJ)) c IF(QN(IJ) .GT. QSIJ) CALL SSHEAT(QN(IJ),RN(IJ),p(ij)) 320 CONTINUE c ---------------------------------------- * PRINTOUT AFTER FIRST TIMESTEP. - * ........................................ * first calculate PBL parameters for - * print out on the file "pblout.dat" - * "work" is the temporary array for - * the potential temperature. - c ---------------------------------------- thbar=0. upbar=0. vpbar=0. do321 i=2,mzh work(i)=rnm(i)*ssk(i) thbar=thbar+work(i) upbar=upbar+unm(i) vpbar=vpbar+vnm(i) 321 continue nlpbl=mzh-1 thbar=thbar/float(nlpbl) upbar=upbar/float(nlpbl) vpbar=vpbar/float(nlpbl) IF(INIT.EQ.0)THEN ihcur=ifix(timeis) ihr=mod(ihcur,24) CALL PRINT INIT=INIT+1 ENDIF c --------------------------------------------------------- * TIMELOOP USING LEAP-FROG TIMESTEPPING FOR MAXSTP STEPS * THEN REINITIALIZING. c --------------------------------------------------------- MAXSTP = 25 DO 380 ITIME=2,MAXSTP delt2=2.*deltat c --------------------------------------------- * compute current value of geostrophic wind. * c --------------------------------------------- tg(1)=timeis ug(1)=xskip vg(1)=xskip call chacha(ug,tg,ns1,ns2,ugi,tgi,mz0,nug,xskip,dz) if(ns1.ne.1.and.ns2.ne.1)call err1(ns1,ns2) call chacha(vg,tg,ns1,ns2,vgi,tgi,mz0,nug,xskip,dz) if(ns1.ne.1.and.ns2.ne.1)call err1(ns1,ns2) c ------------------------------------------------------------ * if precipitation period, let it rain, accumulate rainfall * c ------------------------------------------------------------ IF((TIMEIS.GE.PRST).AND.(TIMEIS.LT.PREND))THEN prcp=prcip else prcp=0. ENDIF ACCP = ACCP + PRCP*DELTAT c ------------------------------------------------------- * solve hydrostatic equation, obtain z at sigma level * c ------------------------------------------------------- CALL SIGTOZ(ZZ,RNM,QNM,MZBL) C ---------------------------------------------------------------------- C ENTER SFC VALUES AND VALUES AT FIRST LEVEL FOR USE IN SFLX. - C ---------------------------------------------------------------------- TS = RNM(1) T2 = RNM(2) cx TH2 = T2 * SS(2)**(-0.286) TH2 = T2 * SSK(2) WS = QNM(1) W2 = QNM(2) U2 = UNM(2) V2 = VNM(2) ZZ2 = ZZ(2) C --------------------------------------- C COMPUTE DOWNWARD RADIATION. C --------------------------------------- TSUN = TIMEIS + TZONE CALL SUN C ---------------------------------------------------------------------- C COMPUTE SURFACE DRAG COEFFICIENTS, UPDATE SOIL, SFC TEMP.,SFC. C MOIST. IF LAST TIMESTEP BEFORE REINITIALIZATION MAKE ONLY 1/2 step C ---------------------------------------------------------------------- IF(ITIME .EQ. MAXSTP)then SOILDT=DELTAT/2. else soildt=deltat ENDIF c CALL SFLX (FDOWN,SSOIL,PSFC,PRCP,DEW,SOILDT,IFHF,IFSNO,ircflg) CALL SFLX . (FDOWN,SSOIL,PSFC,PRCP,DEW,SOILDT,IFHF,IFSNO,ircflg,iflag) C -------------------------------- C ACCUMULATE DEWFALL. C -------------------------------- ACCD = ACCD + DEW*DELTAT C ---------------------------------------------------------------------- C LOAD PBL COLUMN INTO D( )DT ARRAYS, CALL PBL MODEL. ON EXIT C D( )DT ARRAYS CONTAIN TENDENCIES DUE TO PBL MIXING. C ---------------------------------------------------------------------- DO 410 J=2,MZBL DZZ = ZZ(J+1) - ZZ(J) DUDT(J) = UNM(J) DVDT(J) = VNM(J) DTDT(J) = RNM(J)*SSK(J) DWDT(J) = QNM(J) IF (J .LT. MZBL) THEN DUDZ(J) = (UNM(J+1)-UNM(J))/DZZ DVDZ(J) = (VNM(J+1)-VNM(J))/DZZ DQDZ(J) = (QNM(J+1)-QNM(J))/DZZ DTDZ(J) = (RNM(J+1)-RNM(J))/DZZ + G/CP ENDIF 410 CONTINUE CALL PBLL(ifree,envk,ifkpbl,rkfac,rb2) u2t=unm(2) v2t=vnm(2) t2t=rnm(2) thbar=0. upbar=0. vpbar=0. do341 i=2,mzh work(i)=rnm(i)*ssk(i) thbar=thbar+work(i) upbar=upbar+unm(i) vpbar=vpbar+vnm(i) 341 continue c -------------------------------------------------------- * Timestepping: First steps (ITIME=2 and 3) start from * initial fields and do not update previous fields. c * Step 4 also starts from the initial fields, but * updates the previous field. c -------------------------------------------------------- IF(ITIME.le.3)then DO330 IJ=2,MZBL UGZ = UG(1) VGZ = VG(1) uadv=vert(ij)*dudz(ij) vadv=vert(ij)*dvdz(ij) tadv=vert(ij)*dtdz(ij) qadv=vert(ij)*dqdz(ij) fv=fh*(vn(ij)-vgz) fu=fh*(ugz-un(ij)) UP=UNM(IJ)+(fv-uadv+DUDT(IJ))*DELT2 VN(IJ)=VNM(IJ)+(fu-vadv+DVDT(IJ))*DELT2 RN(IJ)=RNM(IJ)+(DTDT(IJ)/SSK(IJ)-tadv)*DELT2 QN(IJ)=QNM(IJ)+(DWDT(IJ)-qadv)*DELT2 c CALL SVP(QSIJ,ESIJ,p(ij),RN(IJ)) QN(IJ) = AMAX1(0.,QN(IJ)) c IF(QN(IJ).GT.QSIJ) CALL SSHEAT(QN(IJ),RN(IJ),p(ij)) UN(IJ) = UP 330 CONTINUE RN(1) = TS QN(1) = WS else DO340 IJ=2,MZBL UGZ = UG(1) VGZ = VG(1) uadv=vert(ij)*dudz(ij) vadv=vert(ij)*dvdz(ij) tadv=vert(ij)*dtdz(ij) qadv=vert(ij)*dqdz(ij) fv=fh*(vn(ij)-vgz) fu=fh*(ugz-un(ij)) UP = UNM(IJ) + (fv - uadv + DUDT(IJ))*DELT2 VP = VNM(IJ) + (fu - vadv + DVDT(IJ))*DELT2 RP = RNM(IJ) + (DTDT(IJ)/SSK(IJ) - tadv)*DELT2 QP = QNM(IJ) + (DWDT(IJ) - qadv)*DELT2 c CALL SVP(QSIJ,ESIJ,p(ij),RP) QP = AMAX1(0.,QP) c IF(QN(IJ).GT.QSIJ) CALL SSHEAT(QN(IJ),RN(IJ),p(ij)) UNM(IJ) = UN(IJ) VNM(IJ) = VN(IJ) RNM(IJ) = RN(IJ) QNM(IJ) = QN(IJ) UN(IJ) = UP VN(IJ) = VP RN(IJ) = RP QN(IJ) = QP 340 CONTINUE rnm(1)=rn(1) qnm(1)=qn(1) rn(1)=ts qn(1)=ws endif c -------------------------------------------------------------------- * Update the model time: * ---------------------- * Multiply time step by 2 for the first two times, * and compute current hour. If the current hour (ihcur) is same as * the one at the previous step (ih2), skip the printing routine. * Reset the latest model time for print out hour. c -------------------------------------------------------------------- timeis=timeis+deltat/3600. if(itime.lt.4)deltat=delt2 ihcur=ifix(timeis) nlpbl=mzh-1 thbar=thbar/float(nlpbl) upbar=upbar/float(nlpbl) vpbar=vpbar/float(nlpbl) if(ihcur.ne.ih2)then ih2=ihcur c ----------------------------------------------------------- * Test for correct modulus of current hour for print out. * (mod of hour interval = 1 here.). If the current model * time is less than the end time, continue the time loop. c ----------------------------------------------------------- if(mod(ihcur,1).eq.0)then ihcur=ifix(timeis) ihr=mod(ihcur,24) call print init=init+1 endif if(timeis.ge.tend)then write (*,*)'# of hours of data written out = ',init close(unit=lunb,status='keep') close(unit=12) stop endif endif 380 CONTINUE c -------------------------------------------------------------------- * Take mean values of n-th and (n-1)-th time step and reinitialize. * This step is made to damp computational mode in leap-frog scheme. * This procedure may be substitited by other mentod such as forward * differencing for every fixed number of leap-frog time steps. c * Note: Soil-canopy model is not stepped backwards since last soil * step was only 1/2 time step. * Atmospheric values after time averaging are referenced to * current time - 1/2 time step c -------------------------------------------------------------------- DO 390 IJ=1,MZBL UN(IJ) = .5 * (UN(IJ)+UNM(IJ)) VN(IJ) = .5 * (VN(IJ)+VNM(IJ)) RN(IJ) = .5 * (RN(IJ)+RNM(IJ)) QN(IJ) = .5 * (QN(IJ)+QNM(IJ)) 390 CONTINUE TIMEIS = TIMEIS - DELTAT/7200. c ----------------------------------------------- * GO TO NEW INITIAL 1/4 FORWARD TIMESTEP. * c ----------------------------------------------- GO TO 301 990 WRITE (*,*) 'END OF FILE REACHED IN CONTROL FILE' WRITE (*,*) ' PREMATURE PROGRAM TERMINATION ' 999 STOP c --------------------------------------------------------------------- * format statements * c --------------------------------------------------------------------- 800 format(A25,$) 810 format(A34) 820 format(/,' CONTROL DATA FILE : ',A34) 830 format(/,' BINARY OUTPUT FILE : ',A34) 840 format(A72) 900 format(1h ,' IFHF, IFCRI, IFSNO, IFCLD =',4i4) 905 format(1h ,' NOOFGR,DELTAT,TEND,RICR,PINK,KOOL =',i4,4f11.3,i4) 910 format(1h ,' TITLE =',a72) 915 format(1h ,' SLA, SLO, TZONE, Z0, Z0H, ZD0, ALBEDO =',7f12.4) 920 format(1h ,' MO, DY, HR =',i3,2f8.2) 925 format(1h ,' PSFC =',f12.2,/,1h ,' TREF =',f12.3,/,1h , . ' CLC =',f12.4,/,1h ,' CMC =',f12.4) 930 format(1h ,'PRST,PREND,PRCIP,ESD,TSNOW =',5f12.4) 935 format(1h ,'ISOIL, TWILT, SIGMAF, IFTC, TSO0, TSOREF, PC =',i3, . 2f12.4,i2,3f12.4) 940 format(1h ,' NUG =',i4) 945 format(1h ,' UG, VG, TG =',3f10.4) 950 format(1h ,' MZ1 =',i4) 955 format(1h ,6f10.3) 960 format(1h ,'NSOIL = ',i3) 965 format(1h ,'DSOIL, WSOIL, TSOIL =',3f12.4) 970 format(1x,'This file written at ',2(i2,'/'),i2,1x,2(i2,':'),i2) 1001 format(1x,'local time and hpbl =',i3,2x,f6.1) 1010 format(1x,f7.4,',',1p,e11.4,',',0p,6(f7.3,','),f6.1) c --------------------------------------------------------------------- END c ------------------------------------------------------------------- c INITIALIZATION OF RADIATION PARAMETERS. c ------------------------------------------------------------------- BLOCK DATA A1 COMMON /AUXI2/ EP,ETOT,PR COMMON /SNOW/ FLX1,FLX2,FLX3 COMMON /RAD/ FDOWN,MO,DY,HO,CLC,TREF,SLO,SLA,ALBEDO,IFCRI,IFPR, . SOLARN COMMON /SFCL/ TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,TSNOW,Z01,CBAGK c ------------------------------------------------------------------- DATA ETOT,PR /0.,1./ DATA FLX1,FLX2,FLX3 /3*0./ DATA TREF,IFPR /270.,1/ DATA TSNOW,CBAGK /273.16,0./ c ---------------------------------------------------------------------- c INITIALIZE OTHER BLOCK DATA AS ZEROS TO SATISFY THE MACINTOSH c MACTRAN+ COMPILER REQUIREMENT THAT ALL COMMON BLOCK DATA BE c INITIALIZED IN THE BLOCK DATA STATEMENTS. THE DATA BELOW WILL BE c RE-INITIALIZED WHEN THE PROGRAM IS EXECUTED AND THE INPUT c (CONTROL) FILE IS READ IN. c ---------------------------------------------------------------------- DATA EP /0./ DATA FDOWN,DY,HO,CLC,SLO,SLA,ALBEDO,SOLARN /8*0./ DATA MO,IFCRI /2*0/ DATA TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,Z01 /17*0./ END c ------------------------------------------------------------------- c INITIALIZATION OF SOIL MODEL PARAMETERS. c ------------------------------------------------------------------- BLOCK DATA A2 PARAMETER (NSOLD=10) LOGICAL KWILT COMMON /SOIL/ WSOIL(NSOLD),CMC,NSOIL,SSOIL,TSOIL(NSOLD),ESD COMMON /SOIL1/ IFTC,TSO0,SIGMAF,SATT,TWILT,SCANOP,PC,RC, . KWILT,CFACTR,TSOREF COMMON /SOIL2/ Z(NSOLD),ISOIL,IFSOIL COMMON /SOIL3/ B(11),SATPSI(11),SATKT(11),TSAT(11) COMMON /SOIL4/ TBOT,ZBOT DATA NSOIL /2/ DATA SATT,KWILT,CFACTR /0.,.FALSE.,0.5/ DATA Z,IFSOIL /-.05,-1.,8*0.,1/ DATA B/4.05,4.38,4.9,5.3,5.39,7.12,7.75,8.52, . 10.4,10.4,11.4/ DATA SATPSI/.121,.09,.218,.786,.478,.299,.356,.63, . .153,.49,.405/ DATA SATKT/1.76E-4,1.5633E-4,3.467E-5,7.2E-6, . 6.95E-6,6.3E-6,1.7E-6,2.45E-6,2.167E-6, . 1.033E-6,1.283E-6/ DATA TSAT/.395,.41,.435,.485,.451,.42,.477,.476, . .426,.492,.482/ DATA TBOT,ZBOT /283.16,-3./ c ---------------------------------------------------------------------- c INITIALIZE OTHER BLOCK DATA AS ZEROS TO SATISFY THE MACINTOSH c MACTRAN+ COMPILER REQUIREMENT THAT ALL COMMON BLOCK DATA BE c INITIALIZED IN THE BLOCK DATA STATEMENTS. THE DATA BELOW WILL BE c RE-INITIALIZED WHEN THE PROGRAM IS EXECUTED AND THE INPUT c (CONTROL) FILE IS READ IN. c ---------------------------------------------------------------------- DATA WSOIL,CMC,SSOIL,TSOIL,ESD /23*0./ DATA TSO0,SIGMAF,TWILT,SCANOP,PC,RC,TSOREF /7*0./ DATA IFTC,ISOIL /2*0/ END c ------------------------------------------------------------------- c INITIALIZATION OF "ALPHA" TABLES FOR SATURATION VAPOR PRESSURE. c ------------------------------------------------------------------- BLOCK DATA A3 COMMON /WFK74/ALF(7,2) DATA ALF/610.7799961,44.36518521,1.428945805, . 2.650648471E-2,3.031240396E-4,2.034080948E-6, . 6.136820929E-9, . 486.6786841,31.52625546,0.8640188586, . 1.279669658E-2,1.077955914E-4,4.886796102E-7, . 9.296950850E-10/ END c ------------------------------------------------------------------- c INITIALIZATION OF BOUNDARY LAYER AND I/O PARAMETERS c ------------------------------------------------------------------- BLOCK DATA A4 PARAMETER (NLEVD=100) COMMON /PBL/ DUDT(NLEVD),DVDT(NLEVD),DTDT(NLEVD),DWDT(NLEVD), . VERT(NLEVD),ZZ(NLEVD),WINIT(NLEVD),HPBL,MZBL,DELT2, . RICR,PINK,DTNW,KOOL,PBLK(NLEVD),CLOUDK(NLEVD), . ZLCL,CLTOP,DCZ,IFCLD COMMON /GRID/ SS(NLEVD),SSK(NLEVD),MZ,MZM,NOOFGR COMMON /XLUN/ LUNC,LUNB,LUND c DATA RICR,PINK,KOOL,PBLK /1.0,2.,0,100*0./ DATA RICR,PINK,KOOL,PBLK /0.25,2.,0,100*0./ DATA NOOFGR /1/ DATA LUNC,LUNB,LUND /2,3,4/ C ---------------------------------------------------------------------- C INITIALIZE OTHER BLOCK DATA AS ZEROS TO SATISFY THE MACINTOSH C MACTRAN+ COMPILER REQUIREMENT THAT ALL COMMON BLOCK DATA BE C INITIALIZED IN THE BLOCK DATA STATEMENTS. THE DATA BELOW WILL BE C RE-INITIALIZED WHEN THE PROGRAM IS EXECUTED AND THE INPUT C (CONTROL) FILE IS READ IN. C ---------------------------------------------------------------------- DATA DUDT,DVDT,DTDT,DWDT,VERT,ZZ,WINIT,CLOUDK /800*0./ DATA SS,SSK,MZ,MZM /200*0.,2*0/ DATA HPBL,DELT2,DTNW,ZLCL,CLTOP,DCZ /6*0./ DATA MZBL,IFCLD /2*0/ END C ********************************************************************** C SUBROUTINE ERR1 C --------------- C WRITES OUT ERROR MESSAGE AND ARRAY BOUNDS PASSED TO "CHACHA" FOR C THE CASE WHERE THE VALUES SHOULD BE 1 BUT ARE CHANGED TO 0 IN C THE SUBROUTINE. C ********************************************************************** SUBROUTINE ERR1(IROW,ICOL) WRITE (*,*)' ERROR IN CALLING SUBROUTINE ''CHACHA''' WRITE (*,*)' THE ROW AND COLUMN DIMENSIONS SHOULD BE=1' WRITE (*,*)' ACTUAL ROW=',IROW,' ACTUAL COLUMN=',ICOL WRITE (*,*)' PROGRAM TERMINATED !!!!!!!!!!' STOP END C ********************************************************************** C SUBROUTINE SFLX C --------------- C COMPUTES DRAG COEFFICIENT FOR MOMENTUM AND HEAT. SECONDLY C POTENTIAL EVAPORATION IS ESTIMATED FROM SURFACE ENERGY BALANCE C (MODIFIED PENMAN). THIRDLY ONE STEP FORWARD IN SOIL MODEL IS C PERFORMED AND ACTUAL EVAPOTRANSPIRATION IS ESTIMATED. THIS C SUBROUTINE WILL ACCOMODATE + OR - SOIL HEAT FLUX. C ---------------------------------------------------------------------- C SOURCE: 5-FEB-87, HUA-LU PAN. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: F, PSFC, PRCP, DT, IFHF, IFSNO C OUTPUT: S, DEW C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SFCL/ T1 (PRESENT), U2, V2, T2, TH2, Q2, Z, Z0, Z0H, T0 C /SOIL/ SMC, CMC, NSOIL, SSOIL, STC C /SOIL2/ ZSOIL, ISOIL C /AUX12/ PR C OUTPUT: C /SFCL/ T1 (UPDATED), Q1, CM, CH, BETA, RIR, RIB, Z01 C /SOIL/ ESD C /AUX12/ EP, ETOT C /SNOW/ FLX1, FLX2, FLX3 C ---------------------------------------------------------------------- C SUBROUTINES CALLED: SFCXCH, HFLX, SVP, SMFLX, KTSOIL, SHFLX, C TDEW C ********************************************************************** SUBROUTINE SFLX(F,S,PSFC,PRCP,DEW,DT,IFHF,IFSNO,ircflg,iflag) c PARAMETER (NSOLD=10) PARAMETER (NSOLD=10,NLEVD=100) COMMON /SFCL/ T1,Q1,CM,CH,U2,V2,T2,TH2,Q2,Z,Z0,Z0H,BETA,RIR, . RIB,RIF,T0,Z01,CBAGK COMMON /SOIL/ SMC(NSOLD),CMC,NSOIL,SSOIL,STC(NSOLD),ESD COMMON /SOIL1/ IFTC,TSO0,SIGMAF,TSOSAT,TWILT,SCANOP, . PC,RC,KWILT,CFACTR,TSOREF INTEGER *4 FLGSN1, FLGSN2 COMMON /SOIL2/ ZSOIL(NSOLD),ISOIL,IFSOIL COMMON /AUXI2/ EP,ETOT,PR COMMON /SNOW/ FLX1, FLX2, FLX3 COMMON /PLANT/ RCMIN COMMON /GRID/ SS(NLEVD),SSK(NLEVD),MZ,MZM,NOOFGR DATA ELV,SIGMA,CP,RGAS,EXMCH /2.5E6,5.67E-8,1004.5,287.,-1./ DATA G /9.806/ DATA DELFAK,SFAK,ELCP /1.3485E7,6.48E-8,2.4888E3/ DATA CPICE /2.106E+3/ DATA ELI,DFS,CPH2O /3.335E5,.13,4.218E3/ C ---------------------------------------------------------------------- C 1. DEFINE COEFFICIENTS. C ---------------------------------------------------------------------- Z01=Z0 S2=U2*U2+V2*V2 SPD=SQRT(S2) T2V=T2*(1.+0.61*Q2) THETAS = T1 * SSK(1) C ---------------------------------------------------------------------- C 2. OVER SEA CONDITIONS ESTIMATE Z0. C ---------------------------------------------------------------------- IF (Z0 .LE. 0.) THEN U3 = S2+1. Z01 = U3*2.1E-6 A = ALOG(Z/Z01) Z01 = 3.E-4*U3/(A*A) ENDIF C ---------------------------------------------------------------------- C 3. COMPUTE BULK RICHARDSON NUMBER AND SURFACE EXCHANGE C COEFFICIENTS. C ---------------------------------------------------------------------- c CALL SFCXCH(T1,T2,TH2,Q1,Q2,Z,Z0,Z0H,Z01,S2,SPD,RIB,CM,CH) CALL SFCXCH . (T1,T2,THETAS,TH2,Q1,Q2,Z,Z0,Z0H,Z01,S2,SPD,RIB,CM,CH) C ---------------------------------------------------------------------- C 4.1 DETERMINE SURFACE MOISTURE (Q1) WHICH DEPENDS ON THE SURFACE C OVER WHICH YOU ARE. C OVER WATER (Z0 = 0) T1 (SEA SURFACE TEMPERATURE OVER THE C OCEAN) IS PRESCRIBED AND Q1 IS THEN THE SATURATION VALUE. C (RETURN TO CALLING PROGRAM AFTER Q1 CALCULATED.) C ---------------------------------------------------------------------- IF (Z0 .EQ. 0) THEN CALL SVP(QSAT,ESAT,PSFC,T1) Q1 = QSAT C ---------------------------------------------------------------------- C 4.2 OVER LAND (Z .NE. 0) USE MODIFIED PENMAN TO ESTIMATE POTENTIAL C EVAPORATION. C ---------------------------------------------------------------------- ELSE TS2=T1*T1 TS4=TS2*TS2 RNET=F-SIGMA*TS4 C ---------------------------------------------------------------------- C 4.2.1 SOIL HEAT FLUX AND NO SNOW. C ---------------------------------------------------------------------- IF ((IFHF .EQ. 1) .AND. (IFSNO .EQ. 0)) THEN CALL HFLX(S,T1,STC,SMC,ZSOIL) C ---------------------------------------------------------------------- C 4.2.2 NO SOIL HEAT FLUX AND NO SNOW. C ---------------------------------------------------------------------- ELSE IF ((IFHF .EQ. 0) .AND. (IFSNO .EQ. 0)) THEN S=0. C ---------------------------------------------------------------------- C 4.2.3 SNOW MODEL: SOIL HEAT FLUX AND SNOW SURFACE (IFSNO=1). C ---------------------------------------------------------------------- ELSE S = 0. FLGSN1 = 0 FLGSN2 = 0 IF (PRCP .GT. 0.) THEN IF (T2 .LE. T0) THEN C ---------------------------------------------------------------------- C 4.2.3.1 FRESH SNOW. C ---------------------------------------------------------------------- FLGSN1=1 ELSEIF (T1.LE.T0) THEN C ---------------------------------------------------------------------- C 4.2.3.2 GROUND IS COLDER THEN 0C, FREEZING RAIN OCCURS. C ---------------------------------------------------------------------- FLGSN2=1 ENDIF ENDIF C ---------------------------------------------------------------------- C INVOKE SNOW MODEL HEAT FLUX. C IF T1 > T0, THERE IS SNOW MELT AND SNOW TEMP WILL BE C T0 (T-ZERO). C ---------------------------------------------------------------------- IF (ESD .GT. 0.) THEN TFLX = AMIN1(T1,T0) S = DFS * (TFLX - STC(1)) / (10. * ESD) C ---------------------------------------------------------------------- C SET MINIMUM FLUX WITHIN SNOW TO -200 W M-2 FOR NUMERICAL C REASONS. C ---------------------------------------------------------------------- S = AMAX1(S,-200.) S = AMIN1(S,200.) ELSE C ---------------------------------------------------------------------- C SOIL HEAT FLUX. C ---------------------------------------------------------------------- CALL HFLX(S,T1,STC,SMC,ZSOIL) ENDIF C ---------------------------------------------------------------------- C CONVERT PRCP TO M/S AND ACCUMULATE IN EQUIVALENT SNOW DEPTH. C ---------------------------------------------------------------------- IF ((FLGSN1 .NE. 0) .OR. (FLGSN2 .NE. 0)) THEN ESD=ESD+PRCP*DT*.001 ENDIF ENDIF T22=T2*T2 T24=T22*T22 CALL SVP(Q2SAT,ESAT,PSFC,T2) DELTA = ELCP * DQSDT(T2,PSFC) RR=T24*SFAK/(PSFC*CH)+1. RHO=PSFC/(RGAS*T2V) RCH=RHO*CP*CH IF (IFSNO .EQ. 1) THEN FLX1=0. C ---------------------------------------------------------------------- C INCLUDE THE WARMING OR COOLING OF RAIN. C ---------------------------------------------------------------------- IF (FLGSN1 .NE. 0) THEN RR = RR + CPICE * PRCP / RCH ELSEIF (PRCP .GT. 0) THEN RR = RR + CPH2O * PRCP / RCH ENDIF ENDIF FNET = F - SIGMA*T24 - S C ---------------------------------------------------------------------- C INCLUDE THE CONVERSION OF RAIN WATER TO ICE. C ---------------------------------------------------------------------- IF (IFSNO .EQ. 1) THEN FLX2=0. IF (FLGSN2.NE.0) THEN FNET=FNET+ELI*PRCP ENDIF IF (FLGSN2.NE.0) THEN FLX2 = -ELI * PRCP ENDIF ENDIF C ---------------------------------------------------------------------- C TERMS USED IN THE CALCULATION OF POTENTIAL EVAPORATION. C EPSCA IS POTENTIAL EVAPORATION IN UNITS OF K M-2 S2; C EP IN KG M-2 S-1. C ---------------------------------------------------------------------- RAD = FNET/RCH + TH2-T2 SPD = AMAX1(SPD,1.E-3) S2 = AMAX1(S2,1.E-3) RIR = G*Z*(-FNET+S)/(T2*SPD*S2) A = ELCP * (Q2SAT-Q2) c EPSCA = (A*RR + RAD*DELTA)/(DELTA+RR) EPSCA = (A*RR + RAD*DELTA + . (A-DELTA*T2)*(SSK(1)-1.0))/(DELTA+RR+SSK(1)-1.0) IF (IFSNO .EQ. 1) THEN FLX3 = 0. ENDIF EP = EPSCA*RCH/ELV C ---------------------------------------------------------------------- C START CALCULATION OF ACTUAL EVAPORATION FOR NO SNOW CASE. C ---------------------------------------------------------------------- IF ((ESD .LE. 0.) .OR. (IFSNO .EQ. 0)) THEN C ---------------------------------------------------------------------- C CHANGE UNIT OF EP TO M/S. C IF EP<0 WE ASSUME DEW TO FORM; DEW ADDED TO PRCP. C CONVERT FROM KG*M-2*S-1 TO M/S. C ---------------------------------------------------------------------- EP1 = EP*0.001 DEW = 0. IF(EP)200,200,210 200 DEW=-EP1 EP1=0. 210 PRCPEF = PRCP*0.001 + DEW C ---------------------------------------------------------------------- C IF ON INPUT PC<0 THEN CANOPY RESISTANCE (RC) IS -PC AND PC IS C THEN CALCULATED FROM RC VIA PENMAN-MONTEITH AT EACH TIME STEP. C OTHERWISE PC IS JUST TAKEN AS SPECIFIED IN THE (INPUT) CONTROL C FILE AND RC IS CALCULATED. C AT THIS TIME THIS CHANGE IS NOT USED FOR THE CASE OF SNOW. C 1E-6 IN THE RC TERM BELOW USED TO AVOID DIVISION BY ZERO FOR C THE CASE WHEN PC=0. C ---------------------------------------------------------------------- IF (PC .LT. 0.) THEN IFLAG = -1 C RC = -PC RCMIN = -PC ENDIF IF (IFLAG .EQ. -1) THEN CALL CANRES (RR,DELTA,ircflg) C PC = (RR + DELTA) / (RR * (1 + RC*CH) + DELTA) ELSE RC = (DELTA + RR) * (1 - PC) / (RR*CH*(PC+1E-6)) ENDIF C ---------------------------------------------------------------------- C DETERMINE THE TOTAL SURFACE MOISTURE FLUX. C ---------------------------------------------------------------------- CALL SMFLX(E1,RUNOFF,SMC,NSOIL,CMC,EP1,DT,PRCPEF) C ---------------------------------------------------------------------- C ACCUMULATE EVAPORATION INTO ETOT. C ---------------------------------------------------------------------- ETOT = ETOT + E1 * DT E = E1*1000. IF(EP)220,230,240 220 BETA=1. E=EP GOTO 250 230 BETA=0. GOTO 250 240 BETA=E/EP 250 CONTINUE C ---------------------------------------------------------------------- C UPDATE SOIL TEMPERATURES & SURFACE TEMPERATURE. C ---------------------------------------------------------------------- IF (IFHF .EQ. 1) THEN CALL KTSOIL(DF1,SMC) YY=T2 + ((F - SIGMA*T24)/RCH+TH2-T2-BETA*EPSCA) / RR ZZ=DF1/(-.5*ZSOIL(1)*RCH*RR) c ZZ1=ZZ+1 ZZ1 = ZZ+(RR-1.0+SSK(1))/RR CALL SHFLX(S,STC,SMC,NSOIL,T1,DT,YY,ZZ1,ZZ) ELSE T1=T2+(RAD-BETA*EPSCA)/RR ENDIF Q1=Q2+E*CP/RCH C ---------------------------------------------------------------------- C ACTUAL EVAPORATION (SUBLIMINATION) OF SNOW; REDUCES SNOW DEPTH. C ---------------------------------------------------------------------- ELSE EP2 = EP*.001*DT BETA = 1. IF (ESD .LT. EP2) THEN BETA = ESD/EP2 ENDIF DEW = 0 IF (EP .LT. 0.) THEN DEW = -EP*.001 ENDIF C ---------------------------------------------------------------------- C RECALCULATE HEAT FLUX FROM/TO ATMOSPHERE TO/FROM SNOW. C ---------------------------------------------------------------------- c T12 = (T2 + ((F-SIGMA*T24)/RCH + TH2-T2 - BETA*EPSCA)/RR c . + DFS*STC(1)/(10.*ESD*RR*RCH))/ c . (1. + DFS/(10.*ESD*RR*RCH)) T12 = (T2 + ((F-SIGMA*T24)/RCH + TH2-T2 - BETA*EPSCA)/RR . + DFS*STC(1)/(10.*ESD*RR*RCH))/ . (1. + (SSK(1)-1.)/RR + DFS/(10.*ESD*RR*RCH)) S = DFS*(T12-STC(1))/(10.*ESD) ESD = AMAX1(0.,ESD-EP2) EP1 = 0. C ---------------------------------------------------------------------- C NO SNOW MELT. C ---------------------------------------------------------------------- IF (T12 .LE. T0) THEN T1 = T12 PRCPEF = 0. C ---------------------------------------------------------------------- C SNOW MELT CONVERTED TO M/S. C ---------------------------------------------------------------------- ELSE EX = .001 * RCH * RR * (T12 - T0) / ELI SNMAX = EX * DT EX = AMIN1(ESD/DT,EX) ESD = AMAX1(ESD - SNMAX, 0.) T1 = AMAX1(T0,T12-ELI*EX*1000./(RCH*RR)) FLX3 = ELI*EX*1000. PRCPEF = EX ENDIF C ---------------------------------------------------------------------- C UPDATE SOIL MODEL, ETOT IN METERS. C ---------------------------------------------------------------------- CALL SMFLX(E1,RUNOFF,SMC,NSOIL,CMC,EP1,DT,PRCPEF) E = BETA * EP ETOT = ETOT + E *.001 *DT ZZ = 0. c ZZ1 = 1. ZZ1 = ZZ+(RR-1.0+SSK(1))/RR YY = STC(1) T11 = T1 CALL SHFLX(S1,STC,SMC,NSOIL,T11,DT,YY,ZZ1,ZZ) Q1 = Q2 + E * CP / RCH IF (FLGSN1 .NE. 0) THEN FLX1 = CPICE*PRCP*(T1-T2) ELSEIF (PRCP .GT. 0) THEN FLX1 = CPH2O*PRCP*(T1-T2) ENDIF ENDIF ENDIF RETURN END C ********************************************************************** C SUBROUTINE LCL C -------------- C THIS ROUTINE TAKES A GIVEN SOUNDING POINT (P,T,Q) AND FINDS THE C LEVEL PLCL WHERE THE MIXING RATIO Q IS THE SATURATION MIXING C RATIO FOR TEMPERATURE TLCL WHERE TLCL/T = (PLCL/P)** RGAS/CP. C ---------------------------------------------------------------------- C SOURCE: CHENG-TSONG CHU 7/8/86, AMENDED BY HUA-LU PAN 6/9/87. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: T, P, Q C OUTPUT: PLCL, TLCL C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C COMMON /WFK74/ ALF(7,2) C ---------------------------------------------------------------------- C SUBROUTINES CALLED: TDEW C ********************************************************************** SUBROUTINE LCL(T,P,Q,PLCL,TLCL) COMMON /WFK74/ ALF(7,2) DATA EPS,RGAS,CP,T0/0.622,287.,1004.5,273.16/ N=0 CALL TDEW(Q,P,T,TG) CPR = CP/RGAS TG = TG - T0 1 IF(TG.LT.-100.) THEN TLCL = 173.15 PLCL = 0. RETURN ENDIF IJ = 1 IF(TG.LT.-50.) IJ = 2 ES = 0. DES = 0. DO 10 K = 7, 1, -1 ES = ALF(K,IJ) + TG * ES IF(K.GT.1) THEN DES= ALF(K,IJ) * (K - 1) + TG * DES ENDIF 10 CONTINUE TG=TG+T0 FT=P*(TG/T)**CPR DFT=(P*CPR*(TG/T)**(CPR-1))/T GT=EPS*ES-Q*FT DGT=EPS*DES-Q*DFT TLCL=TG-GT/DGT N=N+1 IF((Q*1000.).LT.14) THEN IF(ABS(TLCL-TG).LE.0.01) GOTO 2 ELSE IF(ABS(GT*1000./FT).LE.1.0) GOTO 2 ENDIF IF(N.EQ.40) THEN TLCL=AMIN1(TLCL,TG) GOTO 2 ENDIF TG=TLCL-T0 GOTO 1 2 PLCL=P*(TLCL/T)**CPR RETURN END C ********************************************************************** C SUBROUTINE CHACHA C ----------------- C FITS PIECEWISE LINEAR FUNCTION ON ONE GRID TO PIECEWICE LINEAR C FUNCTION ON ANOTHER GRID BY LEAST SQUARES. C C INPUT: C XI,ZI, I= I1,I2 IS FUNCTION VALUES (XI) AND GRID POINT POSITIONS C (ZI) IN THE INTERVAL FROM I1 TO I2. C C INPUT/OUTPUT: C ZO(K),K=MO1,MO2 IS THE GRID ON WHICH THE FUNCTION IS TO BE C ANALYZED ON INPUT. ON OUTPUT MO1,MO2 MAY BE REDEFINED. MO1 AND C MO2 ARE REDEFINED TO THE PART OF THE GRID WHICH OVERLAPS THE C INPUT GRID ZI. IF NO OVERLAP EXISTS THE VALUES MO1=MO2=0 ARE C RETURNED. C C OUTPUT: C XO(K),K=MO1,MO2 THE FUNCTION VALUES ON THE POSSIBLY C REDEFINED GRID. C C IF THE FINAL GRID ZO IS EXTENDING EFFECTIVELY BEYOND ONE OR BOTH C ENDPOINTS ON THE INITIAL GRID THE VALUE OF THE FIRST POINT IN Z0 C OUTSIDE THE ZI GRID IS ADJUSTED TO THE ENDPOINT VALUE IN ZI. IF C MORE THAN ONE POINT FALL OUTSIDE MO1 AND/OR MO2 ARE ADJUSTED. C C IT IS REQUIRED THAT THE TWO GRIDS ZI,ZO CONSISTS OF EVER C INCREASING VALUES OTHERWISE THE RETURN IS MO1=MO2=0. C C THE VARIABLE XSKIP ENABLES THE INCLUSION OF PREDETERMINED C FUNCTION VALUES XO AT ARBITRARY POINTS. ALL POINTS I WHERE XO(I) C ON INPUT IS NOT EQUAL TO XSKIP WILL ON RETURN HAVE THE VALUE C XO(I) PROVIDED ON ENTRY AND POINTS WITH XO(I)=ISKIP WILL BE USED C FOR THE FITTING. THE RETURNED PIECEWISE LINEAR FUNCTION WILL BE C THE LEAST SQUARE FIT WITH THE ADDITIONAL CONSTRAINTS GIVEN BY C PRESET VALUES OF XO(I)'S. C C THE ARRAY Q IS AN INTERNAL WORKARRAY WHICH MUST BE PROVIDED WITH C MO2 ELEMENTS. C ********************************************************************** SUBROUTINE CHACHA(XO,ZO,MO1,MO2,XI,ZI,I1,I2,XSKIP,Q) DIMENSION XO(MO2),ZO(MO2),XI(I2),ZI(I2),Q(MO2) mi1=i1 mi2=i2 IF(MO1.EQ.MO2) GOTO 200 IF(MI1.GE.MI2) GOTO 400 1 IF(ZO(MO1).GE.ZI(MI1)) GOTO 2 MO1=MO1+1 if(mo1.ge.mo2)then mo1=0 mo2=0 return else goto 1 endif c IF(MO1-MO2)1,100,100 2 IF(ZI(MI1+1).GT.ZO(MO1)) GOTO 3 MI1=MI1+1 IF(MI1.EQ.MI2)GOTO 300 GOTO 2 3 IF(MO1.EQ.MO2)GOTO 200 IPI=MI1 IPO=MO1 C=0. CC=C QM=0. UM=0. XM=0. DZ1=1. 10 X=0. Y=0. IPOP=IPO+1 ZL=ZO(IPO) 20 IPIP=IPI+1 IF(ZI(IPIP)-ZO(IPOP))27,28,28 28 ZU=ZO(IPOP) DELZ=ZI(IPIP)-ZI(IPI) if(delz.le.0.)then mo1=0 mo2=0 return endif c IF(DELZ)100,100,21 21 S=(XI(IPIP)-XI(IPI))/DELZ ZU2=ZU*ZU ZU3=ZU2*ZU ZL2=ZL*ZL ZL3=ZL2*ZL F1=S*ZI(IPI)-XI(IPI) F2=ZU2-ZL2 F3=ZU-ZL if(f3.le.0.)then mo1=0 mo2=0 return endif c IF(F3)100,100,22 22 X=X+S/3.*(ZU3-ZL3)+0.5*(XI(IPI)-S*(ZO(IPO)+ZI(IPI)))*F2 . +F1*ZO(IPO)*F3 Y=Y+0.5*S*F2-F1*F3 16 A=C DZ=(ZO(IPOP)-ZO(IPO)) if(dz.le.0.)then mo1=0 mo2=0 return endif c IF(DZ)100,100,31 31 C=DZ/6. CC=C B=2.*(A+C) RHS=XM/DZ1+Y-X/DZ IF(XO(IPO).EQ.XSKIP)GOTO 32 B=1. RHS=XO(IPO) A=0. CC=0. 32 DZ1=DZ XM=X P=A*QM+B QM=-CC/P Q(IPO)=QM UM=(RHS-A*UM)/P XO(IPO)=UM IPO=IPOP IF(ZI(IPIP).GT.ZO(IPOP))GOTO 15 IPI=IPIP IF(IPI.LT.MI2)GOTO 15 MO2=IPOP ZO(MO2)=ZI(MI2) 15 IF(IPO.LT.MO2)GOTO 10 A=C B=2.*A RHS=X/DZ IF(XO(IPO).EQ.XSKIP)GOTO 33 B=1. RHS=XO(IPO) A=0. 33 P=A*QM+B XO(MO2)=(RHS-A*UM)/P DO11 J=MO2-1,MO1,-1 XO(J)=Q(J)*XO(J+1)+XO(J) 11 CONTINUE RETURN 27 ZU=ZI(IPIP) DELZ=ZI(IPIP)-ZI(IPI) if(delz.le.0.)then mo1=0 mo2=0 return endif c IF(DELZ) 100,100,12 12 S=(XI(IPIP)-XI(IPI))/DELZ ZU2=ZU*ZU ZU3=ZU2*ZU ZL2=ZL*ZL ZL3=ZL2*ZL F1=S*ZI(IPI)-XI(IPI) F2=ZU2-ZL2 F3=ZU-ZL if(f3.le.0.)then mo1=0 mo2=0 return endif c IF(F3)100,100,13 13 X=X+S/3.*(ZU3-ZL3)+0.5*(XI(IPI)-S*(ZO(IPO)+ZI(IPI)))*F2 . +F1*ZO(IPO)*F3 Y=Y+0.5*S*F2-F1*F3 ZL=ZU IPI=IPIP IF(IPI.LT.MI2) GOTO 20 MO2=IPOP ZO(MO2)=ZI(MI2) GOTO 16 100 continue MO2=0 MO1=0 RETURN 200 continue IF(MI1.EQ.MI2) GOTO 300 if(zi(mi1).gt.zo(mo1))then mo1=0 mo2=0 return endif c IF(ZI(MI1).GT.ZO(MO1)) GOTO 100 DO201 I=MI1+1,MI2 II=I IF(ZI(I).GE.ZO(MO1)) GOTO 202 201 CONTINUE mo1=0 mo2=0 return 202 continue IF(XO(MO1).NE.XSKIP) RETURN XO(MO1)=(XI(II)-XI(II-1))/(ZI(II)-ZI(II-1))*(ZO(MO1)-ZI(II-1)) . +XI(II-1) RETURN 300 continue if(zo(mo1).ne.zi(mi1))then mo1=0 mo2=0 return endif c IF(ZO(MO1).NE.ZI(MI1)) GOTO 100 XO(MO1)=XI(MI1) RETURN 400 continue if(mi1.gt.mi2)then mo1=0 mo2=0 return endif c IF(MI1.GT.MI2)GOTO 100 if(zo(mo1).gt.zi(mi1))then mo1=0 mo2=0 return endif c IF(ZO(MO1).GT.ZI(MI1))GOTO 100 if(zo(mo2).lt.zi(mi1))then mo1=0 mo2=0 return endif c IF(ZO(MO2).LT.ZI(MI1))GOTO 100 MO2=MO1 ZO(MO1)=ZI(MI1) IF(XO(MO1).NE.XSKIP)RETURN XO(MO1)=XI(MI1) c .................................................. RETURN END C ********************************************************************** C SUBROUTINE SVP C --------------- C THIS ROUTINE CALCULATES SATURATION VAPOR PRESSURE AND SATURATION C SPECIFIC HUMIDITY. C ---------------------------------------------------------------------- C SOURCE: 11 FEBRUARY 1983, HUA-LU PAN, USES A POLYNOMIAL FORMULA C FOLLOWING LOWE AND FICKE (1974). AMENDED 6/9/87 TO INCLUDE TWO C TEMPERATURE RANGES AND TO USE COMMON BLOCK /WFK74/. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: P, TA C OUTPUT: QS, ES C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C COMMON /WFK74/ ALF(7,2) C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE SVP(QS,ES,P,TA) COMMON /WFK74/ ALF(7,2) DATA EPS,T0/.622,273.16/ C ---------------------------------------------------------------------- C THE FORMULA USED CELSIUS TEMPERATURE. C ---------------------------------------------------------------------- T = TA - T0 IF(T.GE.-100.) THEN IJ = 1 IF(T.LE.-50.) IJ = 2 ES = 0. DO 10 K = 7, 1, -1 ES = ES * T + ALF(K,IJ) 10 CONTINUE QS = EPS * ES / P ELSE ES = 0. QS = 0. ENDIF RETURN END C ********************************************************************** C SUBROUTINE ZTOSIG C ----------------- C COMPUTES SIGMA VECTOR CORRESPONDING TO INPUT Z VECTOR. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: T, Q, Z, N C OUTPUT: S C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE ZTOSIG(S,T,Q,Z,N) DIMENSION S(N),T(N),Z(N),Q(N) REAL KAPPA DATA RGAS,CP,G /287.,1004.5,9.806/ KAPPA=RGAS/CP TIV=T(1)*(1.+0.61*Q(1)) S(1)=1. DO 1 I=1,N-1 IP=I+1 TIPV=T(IP)*(1.+0.61*Q(IP)) RAT=TIPV/TIV BK=G*(Z(IP)-Z(I))/(CP*TIV)-0.5+RAT*0.5 GAMMA=-BK+SQRT(BK*BK+RAT) S(IP)=S(I)*GAMMA**(1./KAPPA) TIV=TIPV 1 CONTINUE RETURN END C ********************************************************************** C SUBROUTINE HFAK C --------------- C COMPUTES FACTORS FOR USE IN THE INTEGRATION OF THE HYDROSTATIC C EQUATION: G*(Z(I)-Z(I)) = ALFA(I)*T(I) + BETA(I+1)*T(I+1). C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: S, N C OUTPUT: NONE C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C OUTPUT: C /HYDRO/ ALFA, BETA C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE HFAK(S,N) REAL KAPPA DIMENSION S(N) COMMON /HYDRO/ ALFA(200),BETA(200) DATA RGAS,CP,G /287.,1004.5,9.806/ KAPPA=RGAS/CP DO 1 I=1,N-1 GAMMA=(S(I+1)/S(I))**KAPPA ALFA(I)=0.5*CP*(1.-GAMMA)/G BETA(I+1)=0.5*CP*(1./GAMMA-1.)/G 1 CONTINUE RETURN END C ********************************************************************** C SUBROUTINE SIGTOZ C ----------------- C COMPUTES HEIGHTS AT SIGMA LEVELS FROM HYDROSTATIC EQUATION USING C PRECOMPUTED COEFFICIENTS FROM HFAK. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: T, Q, N C OUTPUT: Z C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /HYDRO/ ALFA, BETA C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE SIGTOZ(Z,T,Q,N) DIMENSION Z(N),T(N),Q(N) COMMON /HYDRO/ ALFA(200),BETA(200) Z(1) = 0. TIMV = T(1) * (1.+0.61*Q(1)) DO 100 I=2,N TIV = T(I) * (1.+0.61*Q(I)) Z(I) = Z(I-1) + (ALFA(I-1)*TIMV + BETA(I)*TIV) TIMV = TIV 100 CONTINUE RETURN END C ********************************************************************** C SUBROUTINE SHFLX C ---------------- C THIS ROUTINE IS USED TO COMPUTE THE SOIL HEAT FLUX. THE SOIL C TEMPERATURE TSOIL IS A DEPENDENT VARIABLE THAT IS UPDATED C THROUGH A PROGNOSTIC EQUATION. THIS VERSION MUST BE USED C TOGETHER WITH A SURFACE ENERGY BALANCE MODEL USING PENMAN C EQUATIONS. C ---------------------------------------------------------------------- C SOURCE: 5-APRIL-1985 HUA-LU PAN C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: TSOIL (CURRENT), THETA, NSOIL, TS (CURRENT), DELTAT, C YY, ZZ1 C OUTPUT: S, TSOIL (UPDATED), TS (UPDATED) C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL2/ Z, ISOIL C /SOIL4/ TBOT, ZBOT C OUTPUT: C /XLUN/ LUNB, LUNC C ---------------------------------------------------------------------- C SUBROUTINES CALLED: HRT, HSTEP, HFLX C ********************************************************************** SUBROUTINE SHFLX(S,TSOIL,THETA,NSOIL,TS,DELTAT,YY,ZZ1,ZZ) PARAMETER (NSOLD=10) DIMENSION TSOIL(NSOIL), THETA(NSOIL), RHSTS(NSOLD) COMMON /SOIL2/ Z(NSOLD), ISOIL, IFSOIL COMMON /SOIL4/ TBOT, ZBOT COMMON /XLUN/ LUNC,LUNB,LUND IF(NSOIL.GT.10) THEN WRITE (*,6001) WRITE (*,*)'PROGRAM TERMINATED !!!!' CLOSE(UNIT=LUNC) CLOSE(UNIT=LUNB,STATUS='KEEP') STOP ENDIF C ---------------------------------------------------------------------- C COMPUTE THE RIGHT HAND SIDE OF DTS/DT EQUATION. C ---------------------------------------------------------------------- CALL HRT(RHSTS,TS,TSOIL,THETA,NSOIL,Z,YY,ZZ1) C ---------------------------------------------------------------------- C UPDATE TSOIL. C ---------------------------------------------------------------------- CALL HSTEP(TSOIL,RHSTS,DELTAT,NSOIL) C ---------------------------------------------------------------------- C COMPUTE SURFACE HEAT FLUX AND SURFACE TEMPERATURE. C ---------------------------------------------------------------------- c TS=(YY+(ZZ1-1.)*TSOIL(1))/ZZ1 TS=(YY+ZZ*TSOIL(1))/ZZ1 CALL HFLX(S,TS,TSOIL(1),THETA(1),Z(1)) 6001 FORMAT(1H1//20X,'***** TOO MANY SOIL LAYERS *****'/ . 20X,'***** MAXIMUM NSOIL IS 10, ABRUPT TERMINATION *****') RETURN END C ********************************************************************** C SUBROUTINE HFLX C --------------- C THIS ROUTINE PROVIDES SOIL HEAT FLUX AT THE SURFACE. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: TS, TS1, T1, Z1 C OUTPUT: S C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: KTSOIL C ********************************************************************** SUBROUTINE HFLX(S,TS,TS1,T1,Z1) CALL KTSOIL(DF,T1) S = DF * (TS1 - TS) / (.5 * Z1 ) RETURN END C ********************************************************************** C SUBROUTINE HRT C -------------- C THIS ROUTINE COMPUTES THE TENDENCY TERMS FOR THE SOIL C THERMODYNAMIC DIFFUSION EQUATION AND THE MATRIX ELEMENTS OF THE C IMPLICIT TIME-INTEGRATION SCHEME. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: TS, TSOIL, THETA, NSOIL, Z, YY, ZZ1 C OUTPUT: RHSTS C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL2/ ISOIL C /SOIL3/ TSAT C /SOIL4/ TBOT, ZBOT C OUTPUT: C /ABCI/ AI, BI, CI C ---------------------------------------------------------------------- C SUBROUTINES CALLED: KTSOIL C ********************************************************************** SUBROUTINE HRT(RHSTS,TS,TSOIL,THETA,NSOIL,Z,YY,ZZ1) PARAMETER (NSOLD=10) DIMENSION RHSTS(NSOIL), TSOIL(NSOIL), THETA(NSOIL), Z(NSOIL) COMMON /SOIL2/ XZ(NSOLD), ISOIL, IFSOIL COMMON /SOIL3/ B(11), SATPSI(11), SATKT(11), TSAT(11) COMMON /SOIL4/ TBOT, ZBOT COMMON /ABCI/ AI(NSOLD), BI(NSOLD), CI(NSOLD) DATA CSOIL/1.26E6/, CH2O/4.2E6/, CAIR/1005./ CALL KTSOIL(DF0,THETA(1)) TZ = THETA(2) IF(THETA(1).GT.THETA(2)) TZ = THETA(1) DTSDZ = (TSOIL(1) - TSOIL(2)) / (-.5 * Z(2)) CALL KTSOIL(DF,TZ) C ---------------------------------------------------------------------- C COMPUTE HEAT CAPACITY OF THE SOIL (INCLUDES CONTRIBUTIONS FROM C WATER, SOIL, AND AIR). C ---------------------------------------------------------------------- HCPCT = THETA(1) * CH2O + (1. - TSAT(ISOIL)) * CSOIL + . (TSAT(ISOIL) - THETA(1)) * CAIR DDZ = 1. / (-.5 * Z(2)) C ---------------------------------------------------------------------- C AI, BI, AND CI ARE THE MATRIX ELEMENTS OF THE TRI-DIAGONAL MATRIX C OF THE IMPLICIT TIME-STEPPING AGORITHM C ---------------------------------------------------------------------- AI(1) = 0. BI(1) = DF * DDZ / (-Z(1) * HCPCT) CI(1) = -BI(1) BI(1) = BI(1) + DF0 / (.5 * Z(1) * Z(1) * HCPCT*ZZ1) S = DF0 * (TSOIL(1) - YY) / (.5 * Z(1)*ZZ1) C ---------------------------------------------------------------------- C RHSTS IS THE TENDENCY TERM. C ---------------------------------------------------------------------- RHSTS(1) = (DF * DTSDZ - S) / (Z(1) * HCPCT) DO 40 K = 2, NSOIL HCPCT = THETA(K) * CH2O + (1. - THETA(K)) * CSOIL IF(K.EQ.NSOIL) GO TO 20 DTSDZ2 = (TSOIL(K) - TSOIL(K+1)) / (.5 * (Z(K-1) - Z(K+1))) TZ2 = THETA(K+1) IF(THETA(K).GT.THETA(K+1)) TZ2 = THETA(K) CALL KTSOIL(DF2,TZ2) DDZ2 = 2. / (Z(K-1) - Z(K+1)) CI(K) = -DF2 * DDZ2 / ((Z(K-1) - Z(K)) * HCPCT) GO TO 30 20 DTSDZ2 = (TSOIL(K) - TBOT) / (.5 * (Z(K-1) + Z(K))-ZBOT) CALL KTSOIL(DF2,THETA(K)) CI(K) = 0. 30 RHSTS(K) = (DF2 * DTSDZ2 - DF * DTSDZ) / ((Z(K) - Z(K-1)) . * HCPCT) AI(K) = -DF * DDZ / ((Z(K-1) - Z(K)) * HCPCT) BI(K) = -(AI(K) + CI(K)) DF = DF2 DTSDZ = DTSDZ2 DDZ = DDZ2 40 CONTINUE RETURN END C ********************************************************************** C SUBROUTINE HSTEP C ---------------- C UPDATE THE SOIL TEMPERATURE FIELD. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: TSOIL (CURRENT), RHSTS (CURRENT), DELTAT, NSOIL C OUTPUT: TSOIL (UPDATED), RHSTS (CURRENT) C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /ABCI/ AI(CURRENT), BI(CURRENT), CI(CURRENT) C OUTPUT: C /ABCI/ AI(UPDATED), BI(UPDATED), CI(UPDATED) C ---------------------------------------------------------------------- C SUBROUTINES CALLED: ROSR12 C ********************************************************************** SUBROUTINE HSTEP(TSOIL,RHSTS,DELTAT,NSOIL) PARAMETER (NSOLD=10) DIMENSION TSOIL(NSOIL), RHSTS(NSOIL) COMMON /ABCI/ AI(NSOLD), BI(NSOLD), CI(NSOLD) DO 10 K = 1, NSOIL RHSTS(K) = RHSTS(K) * DELTAT AI(K) = AI(K) * DELTAT BI(K) = 1. + BI(K) * DELTAT CI(K) = CI(K) * DELTAT 10 CONTINUE CALL ROSR12(CI,AI,BI,CI,RHSTS,RHSTS,NSOIL) DO 20 K = 1, NSOIL TSOIL(K) = TSOIL(K) + CI(K) 20 CONTINUE RETURN END C ********************************************************************** C SUBROUTINE KTSOIL C ----------------- C THIS DETERMINES THE THERMAL CONDUCTIVITY OF THE SOIL. C ---------------------------------------------------------------------- C SOURCE: AL NAKSHABANDI AND KOHNKE (1965). C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: THETA C OUTPUT: DF C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL2/ ISOIL C /SOIL3/ SATPSI,TSAT C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE KTSOIL(DF,THETA) PARAMETER (NSOLD=10) COMMON /SOIL2/ Z(NSOLD), ISOIL, IFSOIL COMMON /SOIL3/ B(11), SATPSI(11), SATKT(11), TSAT(11) DATA KDFKT/0/ C ---------------------------------------------------------------------- C IF KDFKT DOES NOT EQUAL ISOIL THEN THIS IS INITIAL PASS; C EXECUTE THE IF BLOCK, THEN SET KDFKT EQUAL NDKFT SO THIS BLOCK C IS SKIPPED IN SUBSEQUENT PASSES. C ---------------------------------------------------------------------- IF (KDFKT .NE. ISOIL) THEN F1 = ALOG10(SATPSI(ISOIL)) + B(ISOIL) * ALOG10(TSAT(ISOIL)) + 2. KDFKT = ISOIL ENDIF IF (THETA .GT. 0.) THEN PF = F1 - B(ISOIL) * ALOG10(THETA) ELSE PF = 5.2 ENDIF IF (PF .LE. 5.1) DF = EXP(-(2.7 + PF)) * 420. IF (PF .GT. 5.1) DF = .1744 RETURN END C ********************************************************************** C SUBROUTINE SMFLX C ---------------- C THIS ROUTINE IS USED TO COMPUTE THE SOIL MOISTURE FLUX. THE C SOIL MOISTURE CONTENT (PER UNIT VOLUME) THETA IS A DEPENDENT C VARIABLE THAT IS UPDATED AS WELL AS THE CANOPY WATER CONTENT C THROUGH PROGNOSTIC EQUATIONS. C ---------------------------------------------------------------------- C SOURCE: 22 JUL 1983, H.-L. PAN C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: RUNOFF, THETA (CURRENT), NSOIL, CANOPY (CURRENT), C EPOTNL, PRECIP, DELTAT C OUTPUT: EVAP, THETA (UPDATED), CANOPY (UPDATED) C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL1/ IFTC, SIGMAF, TSAT, SCANOP C /SOIL2/ Z, ISOIL, IFSOIL C ---------------------------------------------------------------------- C SUBROUTINES CALLED: THSAT, SEDIR, SET, SEC, SRC, SRT, SSTEP C ********************************************************************** SUBROUTINE SMFLX (EVAP,RUNOFF,THETA,NSOIL,CANOPY,EPOTNL,DELTAT, . PRECIP) LOGICAL KWILT DIMENSION THETA(NSOIL) DIMENSION ET(10), RHSTT(10) COMMON /SOIL1/ IFTC,T0,SIGMAF,TSAT,TWILT,SCANOP,PC,RC,KWILT, . CFACTR,TSOREF COMMON /SOIL2/ Z(10),ISOIL,IFSOIL DATA KDFKT/0/ C ---------------------------------------------------------------------- C CHECK NSOIL. IF NSOIL IS GREATER THAN 10, COMMON /SOIL2/ NEEDS C TO BE CHANGED. C ---------------------------------------------------------------------- IF (NSOIL .GT. 10) THEN WRITE (*,6001) WRITE (*,*) 'PROGRAM TERMINATED !!!!' STOP ENDIF IF (KDFKT .NE. ISOIL) CALL THSAT(TSAT) KDFKT = ISOIL EDIR = 0. EC = 0. ETT = 0. DO 10 K = 1, NSOIL 10 ET(K) = 0. IF (EPOTNL .EQ. 0.) GO TO 30 C ---------------------------------------------------------------------- C SEDIR RETURNS DIRECT EVAPORATION IN EDIR. C ---------------------------------------------------------------------- CALL SEDIR(EDIR,EPOTNL,THETA(1),Z(1)) ETT = 0. EC = 0. C ---------------------------------------------------------------------- C SKIP IF IFTC IS ZERO. C ---------------------------------------------------------------------- IF (IFTC .EQ. 0) GO TO 30 C ---------------------------------------------------------------------- C SET RETURNS TRANSPIRATION IN ET. C ---------------------------------------------------------------------- CALL SET(ET,NSOIL,EPOTNL,THETA,CANOPY,Z) DO 20 K = 1, NSOIL ETT = ETT + ET(K) 20 CONTINUE C ---------------------------------------------------------------------- C SEC RETURNS CANOPY EVAP IN EC. C ---------------------------------------------------------------------- CALL SEC(EC,EPOTNL,CANOPY) 30 CONTINUE EVAP = EDIR + ETT + EC C ---------------------------------------------------------------------- C PCPDRP IS THE COMBINED PRECIP AND DRIP(FROM CANOPY) THAT GOES C INTO THE SOIL. C ---------------------------------------------------------------------- PCPDRP = PRECIP C ---------------------------------------------------------------------- C SRC RETURNS RIGHT HAND SIDE OF CANOPY EQUATION C ---------------------------------------------------------------------- IF (IFTC .EQ. 0) GO TO 40 CALL SRC(RHSCT,EC,PRECIP) C ---------------------------------------------------------------------- C WHEN CANOPY WATER CONTENT EXCEEDS CAPACITY, THE EXCESS IS C DRIPPED TO THE GROUND. C ---------------------------------------------------------------------- DRIP = 0. IF (CANOPY+DELTAT*RHSCT .GT. SCANOP) . DRIP = CANOPY + DELTAT * RHSCT - SCANOP PCPDRP = (1. - SIGMAF) * PRECIP + DRIP / DELTAT 40 CONTINUE C ---------------------------------------------------------------------- C SRT RETURNS RIGHT HAND SIDE OF THETA EQUATION AND ANY RUNOFF C ---------------------------------------------------------------------- CALL SRT(RHSTT,RUNOFF,EDIR,ET,THETA,NSOIL,PCPDRP,Z) C ---------------------------------------------------------------------- C SSTEP RETURNS UPDATED THETA AND CANOPY C ---------------------------------------------------------------------- CALL SSTEP(THETA,CANOPY,RHSTT,RHSCT,DELTAT,NSOIL,RUNOFF) RETURN 6001 FORMAT(1H1//20X,'***** TOO MANY SOIL LAYERS *****'/ . 20X,'***** MAXIMUM NSOIL IS 10, ABRUPT TERMINATION*****') END C ********************************************************************** C SUBROUTINE SEDIR C ---------------- C COMPUTES DIRECT EVAPORATION FROM BARE SOIL. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: EPOTNL, THETA1, Z1 C OUTPUT: EDIR C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL1/ IFTC, TO, SIGMAF C ---------------------------------------------------------------------- C SUBROUTINES CALLED: DFKT C ********************************************************************** SUBROUTINE SEDIR(EDIR,EPOTNL,THETA1,Z1) REAL KT, KT0 LOGICAL KWILT COMMON /SOIL1/IFTC,TSO0,SIGMAF,TSAT,TWILT,SCANOP,PC,RC, . KWILT,CFACTR,TSOREF DATA DF0 /0./ C ---------------------------------------------------------------------- C INITIALIZE DF0 AND KT0 FOR SURFACE MOISTURE EQUAL TO THE AIR DRY C VALUE (TSO0). C ---------------------------------------------------------------------- IF (DF0 .NE. 0.) GO TO 10 CALL DFKT(DF0,KT0,TSO0) 10 CONTINUE C ---------------------------------------------------------------------- C DFKT RETURNS DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY GIVEN THETA. C ---------------------------------------------------------------------- THETA = THETA1 CALL DFKT(DF,KT,THETA) C ---------------------------------------------------------------------- C CALCULATE APPARENT SURFACE SOIL MOISTURE (THSFC) BASED ON C ATMOSPHERIC DEMAND. FX REPRESENTS BARE SOIL EVAPORATION. IF C THSFC IS GREATER THAN OR EQUAL TO THE AIR DRY VALUE (TSO0), THEN C SOIL EVAPORATION IS POTENTIAL (DEMAND CONTROL). OTHERWISE IT IS C CONTROLED BY THE SOIL MOISTURE GRADIENT (FLUX CONTROL); IN THAT C CASE, THSFC IS JUST THE AIR DRY VALUE. C ---------------------------------------------------------------------- THSFC = THETA1 - (-Z1/2)*(EPOTNL + KT)/DF IF (THSFC .GE. TSO0) THEN FX = EPOTNL ELSE FX = DF*(THETA1-TSO0)/(-Z1/2.) - KT FX = AMAX1(FX,0.) THSFC = TSO0 ENDIF C ---------------------------------------------------------------------- C FACTOR IS THE REDUCTION OF DIRECT EVAPORATION DUE TO SHADING BY C VEGETATION. C ---------------------------------------------------------------------- FACTOR = 1. IF (IFTC .NE. 0) FACTOR = (1. - SIGMAF) EDIR = FX * FACTOR RETURN END C ********************************************************************** C SUBROUTINE SEC C -------------- C COMPUTES CANOPY EVAPORATION. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: EPOTNL, CANOPY C OUTPUT: EC C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL1/ SIGMAF, SCANOP, CFACTR C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE SEC(EC,EPOTNL,CANOPY) LOGICAL KWILT COMMON /SOIL1/ IFTC,T0,SIGMAF,TSAT,TWILT,SCANOP,PC,RC, . KWILT,CFACTR,TSOREF EC = SIGMAF * ((CANOPY / SCANOP) ** CFACTR) * EPOTNL RETURN END C ********************************************************************** C SUBROUTINE SRT C -------------- C COMPUTE TIME TENDENCY OF THE SOIL HYDROLOGY DIFFUSION EQUATION C AND THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: EDIR, ET, THETA, NSOIL, PRECIP, Z C OUTPUT: RHSTT, RUNOFF C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL1/ TSAT C OUTPUT: C /ABCI/ AI, BI, CI C ---------------------------------------------------------------------- C SUBROUTINES CALLED: DFKT C ********************************************************************** SUBROUTINE SRT (RHSTT,RUNOFF,EDIR,ET,THETA,NSOIL,PRECIP,Z) DIMENSION RHSTT (NSOIL),ET(NSOIL),THETA(NSOIL),Z(NSOIL) LOGICAL KWILT COMMON /SOIL1/ IFTC,T0,SIGMAF,TSAT,TWILT,SCANOP,PC,RC, . KWILT,CFACTR,TSOREF COMMON /ABCI/ AI(10), BI(10), CI(10) REAL KT,KZ0,INFMAX,KSAT,KT2 DATA DSAT /0./ IF (DSAT .EQ. 0.) CALL DFKT(DSAT,KSAT,TSAT) KZ0 = PRECIP RUNOFF = 0. IF (KZ0 .EQ. 0.) GO TO 10 INFMAX = -DSAT * (TSAT - THETA(1)) / (.5 * Z(1)) + KSAT IF (THETA(1) .GE. TSAT) INFMAX = KSAT IF (KZ0 .LE. INFMAX) GO TO 10 RUNOFF = KZ0 - INFMAX KZ0 = INFMAX 10 TZ = THETA(2) IF (THETA(1) .GT. THETA(2)) TZ = THETA(1) DTDZ = (THETA(1) - THETA(2)) / (-.5 * Z(2)) CALL DFKT(DF,KT,TZ) DDZ = 1. / (-.5 * Z(2)) C ---------------------------------------------------------------------- C AI, BI, AND CI ARE THE TRI-DIAGONAL MATRIX ELEMENTS OF THE C IMPLICIT TIME SCHEME. C ---------------------------------------------------------------------- AI(1) = 0. BI(1) = DF * DDZ / (-Z(1)) CI(1) = -BI(1) C ---------------------------------------------------------------------- C RHSTT IS THE TIME TENDENCY TERM OF THE VOLUMETRIC SOIL WATER C CONTENT EQUATION. C ---------------------------------------------------------------------- RHSTT(1) = (DF * DTDZ + KT - KZ0 + EDIR + ET(1)) / Z(1) DO 40 K = 2, NSOIL IF (K .EQ. NSOIL) GO TO 20 DTDZ2 = (THETA(K) - THETA(K+1)) / (.5 * (Z(K-1) - Z(K+1))) TZ2 = THETA(K+1) IF (THETA(K) .GT. THETA(K+1)) TZ2 = THETA(K) CALL DFKT(DF2,KT2,TZ2) DDZ2 = 2. / (Z(K-1) - Z(K+1)) CI(K) = -DF2 * DDZ2 / (Z(K-1) - Z(K)) GO TO 30 20 DTDZ2 = 0. CALL DFKT(DF2,KT2,THETA(NSOIL)) CI(K) = 0. 30 RHSTT(K) = (DF2*DTDZ2 + KT2 - DF * DTDZ - KT + ET(K)) / . (Z(K) - Z(K-1)) AI(K) = -DF * DDZ / (Z(K-1) - Z(K)) BI(K) = -(AI(K) + CI(K)) IF (K .EQ. NSOIL) GO TO 40 DF = DF2 KT = KT2 DTDZ = DTDZ2 DDZ = DDZ2 40 CONTINUE RETURN END C ********************************************************************** C SUBROUTINE SET C -------------- C COMPUTES TRANSPIRATION DUE TO PLANTS. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: NSOIL, EPOTNL, THETA, CANOPY, Z C OUTPUT: ET C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL1/ SIGMAF, TWILT, SCANOP, PC, KWILT, CFACTR, TSOREF C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE SET(ET,NSOIL,EPOTNL,THETA,CANOPY,Z) DIMENSION ET(NSOIL), THETA(NSOIL), Z(NSOIL) LOGICAL KWILT COMMON /SOIL1/IFTC,T0,SIGMAF,TSAT,TWILT,SCANOP,PC,RC,KWILT, . CFACTR,TSOREF IF (KWILT .OR. THETA(NSOIL) .LT. TWILT) GO TO 20 PART = SIGMAF * PC * EPOTNL * (1. - (CANOPY / SCANOP) ** CFACTR) C ---------------------------------------------------------------------- C GX IS THE SIMPLE TRANSPIRATION FUNCTION. C ---------------------------------------------------------------------- c GX = (THETA(1) - TWILT) / (TSOREF - TWILT) c IF (GX .GT. 1.) GX = 1. c IF (GX .LT. 0.) GX = 0. c ET(1) = (Z(1)/Z(NSOIL)) * GX * PART ET(1) = (Z(1)/Z(NSOIL)) * PART DO 10 K = 2, NSOIL c GX = (THETA(K) - TWILT) / (TSOREF - TWILT) c IF (GX .GT. 1.) GX = 1. c IF (GX .LT. 0.) GX = 0. c ET(K) = ((Z(K)-Z(K-1))/Z(NSOIL)) * GX * PART ET(K) = ((Z(K)-Z(K-1))/Z(NSOIL)) * PART 10 CONTINUE RETURN 20 DO 30 K = 1, NSOIL ET(K) = 0. 30 CONTINUE RETURN END C ********************************************************************** C SUBROUTINE SRC C -------------- C TIME TENDENCY OF THE CANOPY WATER PREDICTION EQUATION IS C COMPUTED HERE. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: EC, PRECIP C OUTPUT: RHSCT C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL1/ SIGMAF C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE SRC (RHSCT,EC,PRECIP) LOGICAL KWILT COMMON /SOIL1/IFTC,T0,SIGMAF,TSAT,TWILT,SCANOP,PC,RC,KWILT, . CFACTR,TSOREF RHSCT = SIGMAF * PRECIP - EC RETURN END C ********************************************************************** C SUBROUTINE SSTEP C ---------------- C UPDATES THE SOIL WATER CONTENT AND THE CANOPY WATER CONTENT. C ---------------------------------------------------------------------- C SOURCE: C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: THETA, CANOPY, RHSTT, RHSCT, DELTAT, NSOIL C OUTPUT: THETA (UPDATED), CANOPY (UPDATED) C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL1/ IFTC, SCANOP C /ABCI/ AI, BI, CI C ---------------------------------------------------------------------- C SUBROUTINES CALLED: ROSR12 C ********************************************************************** SUBROUTINE SSTEP(THETA,CANOPY,RHSTT,RHSCT,DELTAT,NSOIL,RUNOFF) DIMENSION THETA(NSOIL), RHSTT(NSOIL) LOGICAL KWILT COMMON /SOIL1/IFTC,T0,SIGMAF,TSAT,TWILT,SCANOP,PC,RC,KWILT, . CFACTR,TSOREF COMMON /ABCI/AI(10), BI(10), CI(10) DO 10 K = 1, NSOIL RHSTT(K) = RHSTT(K) * DELTAT AI(K) = AI(K) * DELTAT BI(K) = 1. + BI(K) * DELTAT CI(K) = CI(K) * DELTAT 10 CONTINUE CALL ROSR12(CI,AI,BI,CI,RHSTT,RHSTT,NSOIL) DO 20 K = 1, NSOIL THETA(K) = THETA(K) + CI(K) IF (THETA(K) .LT. 0.) THEN THETA(K) = 0. ENDIF IF (THETA(K) .GT. TSAT) THEN c RUNOFF = RUNOFF + (THETA(K) - TSAT) THETA(K) = TSAT ENDIF 20 CONTINUE IF (IFTC .EQ. 0) RETURN CANOPY = CANOPY + DELTAT * RHSCT IF (CANOPY .LT. 0.) CANOPY = 0. IF (CANOPY .GT. SCANOP) CANOPY = SCANOP RETURN END C ********************************************************************** C SUBROUTINE DFKT C --------------- C COMPUTE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR C THE 11 USDA SOIL TYPES (SEE USER'S GUIDE FOR LIST). C ---------------------------------------------------------------------- C SOURCE: CLAPP AND HORNBERGER C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: THETA C OUTPUT: DF, KT C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL2/ ISOIL, IFSOIL C /SOIL3/ B , SATPSI , SATKT , TSAT C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE DFKT(DF,KT,THETA) REAL KT,KTK COMMON /SOIL2/ Z(10),ISOIL,IFSOIL COMMON /SOIL3/ B(11),SATPSI(11),SATKT(11),TSAT(11) DIMENSION DFK(22),KTK(22) DATA KDFKT /0/ IF (KDFKT .EQ. ISOIL) GO TO 20 F1 = B(ISOIL) * SATKT(ISOIL) * SATPSI(ISOIL) / . TSAT(ISOIL) ** (B(ISOIL) + 3.) F2 = SATKT(ISOIL) / TSAT(ISOIL) ** (B(ISOIL) * 2. + 3.) KDFKT = ISOIL IF (IFSOIL .EQ. 1) GO TO 20 DYNW = TSAT(ISOIL) *.05 DO 10 K = 1, 22 DFK(K) = F1 * (FLOAT(K-1) * DYNW) ** (B(ISOIL) + 2.) KTK(K) = F2 * (FLOAT(K-1) * DYNW) ** (B(ISOIL) * 2. + 3.) 10 CONTINUE 20 CONTINUE IF (IFSOIL .EQ. 2) GO TO 30 DF = F1 * THETA ** (B(ISOIL) + 2.) KT = F2 * THETA ** (B(ISOIL) * 2. + 3.) RETURN 30 W = THETA / DYNW KW = W DF = DFK(KW) + (W - KW) * (DFK(KW+1) - DFK(KW)) KT = KTK(KW) + (W - KW) * (KTK(KW+1) - KTK(KW)) RETURN END C ********************************************************************** C SUBROUTINE THSAT C ---------------- C DETERMINES SOIL MOISTURE SATURATION VALUE FOR SOIL TYPE ISOIL C ---------------------------------------------------------------------- C SOURCE: STANDARD METEOROLOGICAL PRINCIPLES. C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: NONE C OUTPUT: THETA C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /SOIL2/ ISOIL C /SOIL3/ TSAT C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE THSAT(THETA) COMMON /SOIL2/ Z(10),ISOIL,IFSOIL COMMON /SOIL3/ B(11),SATPSI(11),SATKT(11),TSAT(11) THETA = TSAT(ISOIL) RETURN END C ********************************************************************** C SUBROUTINE TDEW C --------------- C DEW-POINT TEMPERATURE IS OBTAINED VIA ADIABATIC-ISOBARIC COOLING C PROCEDURE UNTIL THE MIXING RATIO BECOMES THE SATURATION MIXING C RATIO. C -- C MODIFICATIONS: NONE. C ---------------------------------------------------------------------- C SOURCE: STANDARD METEOROLOGICAL PRINCIPLES C ---------------------------------------------------------------------- C SUBROUTINE ARGUMENTS: C INPUT: Q, P, TA C OUTPUT: TD, C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C /WFK74/: ALF(7,2) C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE TDEW(Q,P,TA,TD) COMMON /WFK74/ ALF(7,2) DATA EPS,T0 /.622,273.16/ TN = TA - T0 IF (Q .LE. 0.) THEN TD = T0 - 99.9 RETURN ENDIF EVP = P * Q/EPS 1 IF (TN .LT. -80.) THEN C ---------------------------------------------------------------------- C CONVERGENCE IS UNCERTAIN IF TN < -80 C ---------------------------------------------------------------------- TD = T0 - 99.9 RETURN ENDIF IJ = 1 IF (TN .LT. -50.) IJ = 2 ES = 0. DFN = 0. DO 10 K = 7, 1, -1 ES = ES * TN + ALF(K,IJ) IF (K .GT. 1) THEN DFN = ALF(K,IJ) * (K - 1) + TN * DFN ENDIF 10 CONTINUE FN = ES - EVP TN1 = TN - FN/DFN IF (ABS(TN1-TN) .LE. 0.01) GOTO 2 TN = TN1 GOTO 1 2 TD = TN1 + T0 RETURN END C ********************************************************************** C SUBROUTINE SP C ------------- C DESCRIPTION: COMPUTES GEOPOTENTIAL HEIGHT OF THE LCL LEVEL. C ---------------------------------------------------------------------- C SOURCE: STANDARD METEOROLOGICAL PRINCIPLES. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: SS, PSFC, ZZ, MZ, NLEVD C OUTPUT: PLCL, ZSP C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE SP(SS,PSFC,PLCL,ZZ,J,MZ,ZSP,NLEVD) DIMENSION SS(NLEVD),ZZ(NLEVD) C ---------------------------------------------------------------------- C RETURN ZSP = -1 IF HEIGHT OF LCL IS AT THE SURFACE. C ---------------------------------------------------------------------- ZSP = -1. IF (PLCL .LE. (PSFC*SS(MZ))) RETURN IF (PLCL .GE. PSFC) THEN ZSP = 0. RETURN ENDIF DO 10 I = MZ-1, 1 ,-1 P1 = PSFC * SS(I) IF (PLCL .LE. P1) GOTO 1 10 CONTINUE 1 P2 = PSFC * SS(I+1) RATIO = (P1-PLCL)/(P1-P2) ZSP=(ZZ(I+1)-ZZ(I))*RATIO+ZZ(I) RETURN END C ********************************************************************** C FUNCTION DQSDT C -------------- C THIS FUNCTION CALCULATES DQSDT = EPS * DESDT / P, TO BE USED IN C NEWTON-RAPHSON ITERATION INVOLVING QS. C ---------------------------------------------------------------------- C SOURCE: STANDARD METEOROLOGICAL PRINCIPLES. C ---------------------------------------------------------------------- C FUNCTION ARGUMENTS: C INPUT: TA, P C OUTPUT: DSQDT C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS CALLED: C INPUT: C /WFK74/ ALF(7,2) C ---------------------------------------------------------------------- C SUBROUTINES USED: NONE C ********************************************************************** FUNCTION DQSDT(TA,P) COMMON /WFK74/ ALF(7,2) DATA EPS,T0 /.622,273.16/ C ---------------------------------------------------------------------- C CONVERSION TO T CELSIUS. C ---------------------------------------------------------------------- T = TA - T0 IF(T.GE.-100.) THEN IJ = 1 C ---------------------------------------------------------------------- C IF TEMPERATURE LESS THAN -50 C, USE SECOND SET OF ALF TABLES. C ---------------------------------------------------------------------- IF (T. LT. -50.) IJ = 2 DESDT = 0. C ---------------------------------------------------------------------- C COMPUTE DES/DT. C ---------------------------------------------------------------------- DO 10 K = 7, 2, -1 DESDT = ALF(K,IJ) * (K - 1) + T * DESDT 10 CONTINUE C ---------------------------------------------------------------------- C CALCULATE DQSDT. C ---------------------------------------------------------------------- DQSDT = EPS * DESDT / P ELSE C ---------------------------------------------------------------------- C IF TEMPERATURE LESS THAN -100 C, SET DQSDT EQUAL TO ZERO. C ---------------------------------------------------------------------- DQSDT = 0. ENDIF RETURN END C ********************************************************************** C SUBROUTINE TOPQ C --------------- C THIS ROUTINE PROVIDES T2 AND QSAT AT T2 AT PRESSURE P2 THAT C GIVES THE SAME EQUIVALENT POTENTIAL TEMPERATURE AS THE POINT C (T1, P1). C ---------------------------------------------------------------------- C SOURCE: STANDARD METEOROLOGICAL PRINCIPLES C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: P2, Q1, T1, P1 C OUTPUT: T2, Q2 C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES USED: SVP C ********************************************************************** SUBROUTINE TOPQ(Q2,T2,P2,Q1,T1,P1) REAL LV,KAPPA DATA CP,DTDP,G,LV,RGAS /1004.5,4.5E-4,9.806,2.5E6,287./ FUNC(QS,T) = EXP(LV * QS / (CP * T)) KAPPA = RGAS / CP C ---------------------------------------------------------------------- C FIRST GUESS. C ---------------------------------------------------------------------- T2 = T1 + DTDP * (P2 - P1) PFACT = (1.E5 / P2) ** KAPPA CONST = T1 * (1.E5 / P1) ** KAPPA * FUNC(Q1,T1) C ---------------------------------------------------------------------- C ITERATION STARTS. C ---------------------------------------------------------------------- ITER = 0 10 CALL SVP(Q2,E2,P2,T2) FACT1 = FUNC(Q2,T2) F = T2 * PFACT * FACT1 - CONST DFDT = PFACT * FACT1 + PFACT * FACT1 * (LV * DQSDT(T2,P2) / . CP - LV * Q2 / (CP * T2)) DT = -F / DFDT T2 = T2 + DT C ---------------------------------------------------------------------- C STOP ITERATIONS IF DT LESS THAN 0.1, OTHERWISE CONTINUE. C ---------------------------------------------------------------------- IF (ABS(DT) .LT. .1) GOTO 100 ITER = ITER + 1 IF (ITER .LT. 100) GOTO 10 STOP 'NO CONVERGENCE IN TOPQ' 100 CALL SVP(Q2,E2,P2,T2) RETURN END C ********************************************************************** C SUBROUTINE SSHEAT C ----------------- C THIS ROUTINE CALCULATES THE SUPER-SATURATION ADJUSTMENT. C ---------------------------------------------------------------------- C SOURCE: STANDARD METEOROLOGICAL PRINCIPLES C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: Q (SUPERSATURATED, Q GREATER THAN QSAT(T)), T, P C OUTPUT: Q (ADJUSTED), T (ADJUSTED) C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: SVP C ********************************************************************** SUBROUTINE SSHEAT (Q,T,P) REAL LV DATA CP,LV /1004.5,2.5E6/ C ---------------------------------------------------------------------- C FIRST GUESS. C ---------------------------------------------------------------------- T2 = T C ---------------------------------------------------------------------- C START ITERATION. C ---------------------------------------------------------------------- ITER = 0 CALL SVP(Q2,E2,P,T2) IF (Q2 .GE. Q) RETURN C ---------------------------------------------------------------------- C CONTINUE ITERATION. IF ERROR CRITERIA SATISFIED, GO TO 20 AND C CALL SVP TO FIND ADJUSTED Q AND T. C ---------------------------------------------------------------------- 10 F = T2 - T - LV * (Q - Q2) / CP DFDT = 1. + LV * DQSDT(T2,P) / CP DT = - F / DFDT T2 = T2 + DT IF (ABS(DT) .LE. .01) GOTO 20 ITER = ITER + 1 CALL SVP(Q2,E2,P,T2) C ---------------------------------------------------------------------- C CONTINUE ITERATIONS IF THE NUMBER OF ITERATIONS IS LESS THAN C 100. C ---------------------------------------------------------------------- IF (ITER .LT. 100) GOTO 10 STOP ' ITERATION IN SSHEAT DID NOT CONVERGE' 20 CALL SVP(Q,E2,P,T2) T = T2 RETURN END C ********************************************************************** C SUBROUTINE SFCXCH C ----------------- C THIS ROUTINE CALCULATES THE BULK RICHARDSON NUMBER AND THE C SURFACE EXCHANGE COEFFICIENTS FOR MOMENTUM AND HEAT (=MOISTURE). C ---------------------------------------------------------------------- C SOURCES: LOUIS, ET AL., 1982; HOLTSLAG AND BELJAARS, 1989. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: T1, T2, TH2, Q2, Z, Z0, Z0H, Z01, S2, SPD C OUTPUT: RIB, CM, CH C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** c SUBROUTINE SFCXCH(T1,T2,TH2,Q1,Q2,Z,Z0,Z0H,Z01,S2,SPD,RIB,CM,CH) SUBROUTINE SFCXCH . (T1,T2,THETAS,TH2,Q1,Q2,Z,Z0,Z0H,Z01,S2,SPD,RIB,CM,CH) PARAMETER (NSOLD=10) DATA PRNEUT,EXMCH /1.,-1./ C DATA B1,B2,B3,B4,CUS,VK,G /10.,15.,10.,8.,7.5,0.40,9.806/ C DATA B1,B2,CUS,RCUCT,VK,G /9.4,9.4,7.4,0.716,0.40,9.806/ C DATA B1,B2,CUS,RCUCT,VK,G /10.,15.,7.5,0.716,0.40,9.806/ DATA BUS,BTS,CUS,CTS,VK,G /10.,15.,7.5,5.,0.40,9.806/ C ---------------------------------------------------------------------- C INITIALIZE Z01, TSV, AND TH2V. C ---------------------------------------------------------------------- Z01=Z0 c TSV=T1*(1.+0.61*Q1) c THSV=THETAS*(1.+0.61*Q1) THSV=THETAS*(1.+0.61*Q2) TH2V=TH2*(1.+0.61*Q2) C ---------------------------------------------------------------------- C COMPUTE BULK RICHARDSON NUMBER AND CONSTANTS. C ---------------------------------------------------------------------- AA = VK/ALOG(Z/Z01) AH = AA*VK/ALOG(Z/Z0H) AM = AA*AA IF (SPD .NE. 0.) THEN c RIB = G*Z*(TH2V-TSV)/(TH2V*S2) RIB = G*Z*(TH2V-THSV)/(TH2V*S2) IF(RIB)140,130,130 C ---------------------------------------------------------------------- C STABLE STRATIFICATION USING METHOD OF MAHRT (1987). C IF RIB IS SO LARGE THAT THE NUMERICAL CALCULATION OF CM = 0, C SPECIFY CM AND CH AS A VERY SMALL NUMBER. C ---------------------------------------------------------------------- 130 CM = AM*SPD*EXP(EXMCH*RIB) IF (CM .LT. 1.E-6) THEN CM = 1.E-6 CH = 1.E-6 ELSE CH = (AH*SPD*EXP(EXMCH*RIB))/PRNEUT ENDIF GOTO 150 C ---------------------------------------------------------------------- C UNSTABLE STRATIFICATION FROM LOUIS, ET AL. (1982). C ---------------------------------------------------------------------- 140 CUP = CUS*AM*BUS*SQRT(-RIB*Z/Z01) CTP = CTS*AH*BTS*SQRT(-RIB*Z/Z0H) CM = (1. - BUS*RIB/(1.+CUP)) * AM*SPD CH = (1. - BTS*RIB/(1.+CTP)) * AH*SPD 150 CONTINUE ELSE C ---------------------------------------------------------------------- C NO WIND; DETERMINE IF UNSTABLE OR STABLE CASE. C ---------------------------------------------------------------------- c DTV = TH2V - TSV c IF (DTV .LT. 0.) THEN DTHV = TH2V - THSV IF (DTHV .LT. 0.) THEN C ---------------------------------------------------------------------- C NO WIND, UNSTABLE (FREE CONVECTION) CASE. C ---------------------------------------------------------------------- c CM = SQRT(-G*DTV*Z01/TSV) / CUS c CH = SQRT(-G*DTV*Z0H/TSV) / CTS CM = SQRT(-G*DTHV*Z01/THSV) / CUS CH = SQRT(-G*DTHV*Z0H/THSV) / CTS ELSE C ---------------------------------------------------------------------- C NO WIND, STABLE CASE; LET RIB = A VERY LARGE NUMBER. (THE C VALUES BELOW ARE NECESSARY TO AVOID LATER DIVISION BY ZERO.) C ---------------------------------------------------------------------- CM = 1.E-6 CH = 1.E-6 RIB = 1.E3 ENDIF ENDIF RETURN END C ********************************************************************** C FUNCTION T2M C ------------ C THIS FUNCTION CALCULATES THE POTENTIAL TEMPERATURE AT 2 METERS. C ---------------------------------------------------------------------- C SOURCE: HOLTSLAG,1987,KNMI SCI REP 87-2. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: PZ, PT2, PHABS, PTST, PUST, PVIRFX C OUTPUT: T2M C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** FUNCTION T2M(PZ,PT2,PHABS,PTST,PUST,PVIRFX) REAL K,G,OBK PARAMETER(K=0.40,G=9.81,OBK=99999.) IF (PUST.NE.OBK.AND.PTST.NE.OBK.AND.PHABS.NE.OBK) THEN IF (PHABS.LT.100) THEN PHABS=PHABS+273.15 ENDIF IF (PVIRFX.EQ.0.0) THEN C ---------------------------------------------------------------------- C NEUTRAL C ---------------------------------------------------------------------- OBUK2=-1E6 ELSE C ---------------------------------------------------------------------- C STABLE OR UNSTABLE C ---------------------------------------------------------------------- OBUK2= -(PHABS*PUST**3)/(PVIRFX*G*K) ENDIF IF ((PTST.EQ.0.).AND.(PUST.LE.0.01)) THEN C ---------------------------------------------------------------------- C VERY STABLE C ---------------------------------------------------------------------- OBUK2=1E-6 ENDIF ENDIF T2M=PT2-(PTST/K*(ALOG(PZ/2)-FPSIH(PZ/OBUK2)+FPSIH(2/OBUK2))) RETURN END C ********************************************************************** C REAL FUNCTION FPSIM C ------------------- C THIS FUNCTION CALCULATES THE STABILITY CORRECTION FUNCTION IN C THE SURFACE LAYER WIND PROFILE. C ---------------------------------------------------------------------- C SOURCE: THE PRESENT MODEL IS AN EMPIRICAL FIT BY HOLTSLAG AND C DE BRUIN(1988, JAM, 22,P689-704) OF DATA BY HICKS(1976, QUART. C J. R. METEOR. SOC., 102, 535-551); SEE ALSO HOLTSLAG(1984, BLM, C 29, 225-250). C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: ETA C OUTPUT: FPSIM C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** REAL FUNCTION FPSIM(ETA) IMPLICIT CHARACTER*1 (A-Z) REAL ETA REAL X, PID2 DATA PID2/1.5707963/ IF (ETA .LT. 0.) THEN X=SQRT(SQRT(1-16.*ETA)) FPSIM=ALOG((1+X)**2*(1+X**2)/8) -2.*ATAN(X) +PID2 ELSE IF (ETA .GT. 200) THEN FPSIM=-0.7*ETA -10.72 ELSE FPSIM=-0.7*ETA -(0.75*ETA-10.72)*EXP(-0.35*ETA) -10.72 ENDIF ENDIF RETURN END C ********************************************************************** C REAL FUNCTION FPSIH C ------------------- C THIS FUNCTION CALCULATES THE STABILITY CORRECTION FUNCTION IN C THE SURFACE LAYER TEMPERATURE PROFILE. C ---------------------------------------------------------------------- C SOURCE: THE PRESENT MODEL IS AN EMPIRICAL FIT BY HOLTSLAG AND C DE BRUIN(1988, JAM,22,P689-704) OF DATA BY HICKS(1976, QUART. C J. R. METEOR. SOC., 102, 535-551); SEE ALSO HOLTSLAG(1984, BLM, C 29, 225-250). C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: ETA C OUTPUT: FPSIH C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: NONE C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** REAL FUNCTION FPSIH(ETA) IMPLICIT CHARACTER*1 (A-Z) REAL ETA REAL Y IF (ETA .LT. 0.) THEN Y=SQRT(1-16.*ETA) FPSIH=2*ALOG((1+Y)/2) ELSE FPSIH=-0.7*ETA -(0.75*ETA-10.72)*EXP(-0.35*ETA) -10.72 ENDIF RETURN END C ---------------------------------------------------------------------- C SUBROUTINE ROSR12 C ----------------- C ROSR12 inverts the tri-diagonal matrix shown in the below: C C ### ### ### ### ### ### C #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # # C #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # # C # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) # C # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) # C # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) # C # . . # # . # = # . # C # . . # # . # # . # C # . . # # . # # . # C # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)# C # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)# C # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) # C ### ### ### ### ### ### c----------------------------------------------------------------------- c DELTA IS A WORKING ARRAY OF DIMENSION M. IF THE ARRAY D IS NOT c REQUIRED AFTERWARD, THE CALL: CALL ROSR12(P,A,B,C,D,D,M), c WILL REDUCE STORAGE REQUIREMENTS. c c IF THE ARRAY C IS NOT REQUIRED AFTERWARD, THE CALL: c CALL ROSR12(C,A,B,C,D,D,M), WILL FURTHER REDUCE STORAGE c REQUIREMENTS. c-------------------- c Modified: 1 March 1990 by Jinwon Kim c * Internal calculation is performed by double-precision c variables. c * Internal array has been increased to 100 (can handle c models upto 100 model layers) c ----------------------------------------------------------------- c ARGUMENTS: c INPUT: aa, bb, cc, dd, m c OUTPUT: pp, delt c ----------------------------------------------------------------- subroutine rosr12 (pp,aa,bb,cc,dd,delt,m) parameter(nl=100) dimension pp(nl),aa(nl),bb(nl),cc(nl),dd(m),delt(m) double precision p(nl),a(nl),b(nl),c(nl),d(nl),delta(nl) double precision x do1 i=1,m a(i)=aa(i) b(i)=bb(i) c(i)=cc(i) d(i)=dd(i) 1 continue c(m)=0. x=1./b(1) p(1)=-c(1)*x delta(1)=d(1)*x do10 k=2,m km=k-1 x=1./(b(k)+a(k)*p(km)) p(k)=-c(k)*x delta(k)=(d(k)-a(k)*delta(km))*x 10 continue p(m)=delta(m) mp=m+1 do20 k=2,m kk=mp-k p(kk)=p(kk)*p(kk+1)+delta(kk) 20 continue do30 i=1,m pp(i)=p(i) 30 continue return end C ---------------------------------------------------------------------- C SUBROUTINE PRINT C ---------------- C ALL BINARY OUTPUT IS WRITTEN FROM THIS ROUTINE. C ---------------------------------------------------------------------- SUBROUTINE PRINT PARAMETER (NLEVD=100,NLEVI=100,NSOLD=10) COMMON /FIELDS/ UN(NLEVD),UNM(NLEVI),VN(NLEVD),VNM(NLEVI), . RN(NLEVD),RNM(NLEVI),QN(NLEVD),QNM(NLEVI), . UG(NLEVD),VG(NLEVD) COMMON /GRID/ SS(NLEVD),SSK(NLEVD),MZ,MZM,NOOFGR COMMON /SOIL/ WSOIL(NSOLD),CMC,NSOIL,SSOIL,TSOIL(NSOLD),ESD COMMON /SOIL1/ IFTC,TSO0,SIGMAF,TSOSAT,TWILT,SCANOP, . PC,RC,KWILT,CFACTR,TSOREF COMMON /SOIL2/ DSOIL(NSOLD),ISOIL,IFSOIL COMMON /PBL/ DUDT(NLEVD),DVDT(NLEVD),DTDT(NLEVD),DWDT(NLEVD), . VERT(NLEVD),ZZ(NLEVD),WINIT(NLEVD),HPBL,MZBL,DELT2, . RICR,PINK,DTNW,KOOL,PBLK(NLEVD),CLOUDK(NLEVD), . ZLCL,CLTOP,DCZ,IFCLD COMMON /SFCL/ TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,TSNOW,Z01,CBAGK COMMON /AUXIL/ PSFC,HEAT,TAOX,TAOY,EVAP,XL,CG,THERM,DELTAT,FH, . PRCP,DEW,ACCP,ACCD,UFLX(NLEVD),VFLX(NLEVD), . TFLX(NLEVD),QFLX(NLEVD) COMMON /RAD/ FDOWN,MO,DY,TSUN,CLC,TREF,SLO,SLA,ALBEDO,IFCRI,IFPR, . SOLARN COMMON /AUXI2/ EP,ETOT,PR COMMON /SNOW/ FLX1, FLX2, FLX3 COMMON /XLUN/ LUNC,LUNB,LUND COMMON /HCALC/ JPBL,HO dimension hfl(nlevd),vfl(nlevd),ufl(nlevd),wfl(nlevd),pwfl(nlevd) dimension zsp(nlevd),tcl(nlevd),plcl(nlevd),cldlqd(nlevd) dimension ptnlqd(nlevd),rh(nlevd),qsat(nlevd) c ----------------------------------------------------------------------- one3d=1./3. ido=ho/24. xho=ho-float(24*ido) iho=xho xmin=60.*(xho-float(iho)) min=xmin sec=60*(xmin-float(min)) write(lunb)mo,ifix(dy),ido,iho,min,sec WRITE (*,*) . '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>' write (*,*) 'DAY,HOUR=',IDO,IHO c ----------------------- c SURFACE FLUXES. c ----------------------- UFL(1) = TAOX VFL(1) = TAOY HFL(1) = HEAT WFL(1) = EVAP c --------------------------------------- c CONSTANT FLUX LAYER IS ASSUMED. c --------------------------------------- UFL(2) = UFL(1) VFL(2) = VFL(1) HFL(2) = HFL(1) WFL(2) = WFL(1) DO 100 J=2, MZM DZ=ZZ(J+1)-ZZ(J) JP=J+1 UFL(JP)=UFL(J)-DZ*0.5*(UFLX(JP)+UFLX(J)) VFL(JP)=VFL(J)-DZ*0.5*(VFLX(JP)+VFLX(J)) HFL(JP)=HFL(J)-DZ*0.5*(TFLX(JP)+TFLX(J)) WFL(JP)=WFL(J)-DZ*0.5*(QFLX(JP)+QFLX(J)) 100 CONTINUE UU = UFL(MZBL) VV = VFL(MZBL) HH = HFL(MZBL) WW = WFL(MZBL) DO 110 J = MZBL, 3, -1 DU = UU - UFL(J) DV = VV - VFL(J) DH = HH - HFL(J) DW = WW - WFL(J) IF(DU.EQ.0.) UFL(J) = 0. IF(DV.EQ.0.) VFL(J) = 0. IF(DH.EQ.0.) HFL(J) = 0. IF(DW.EQ.0.) WFL(J) = 0. 110 CONTINUE PPBL = PSFC*(SS(JPBL)+(HPBL-ZZ(JPBL))*(SS(JPBL+1)-SS(JPBL))/ . (ZZ(JPBL+1)-ZZ(JPBL))) DO 200 J = 1, MZBL P = PSFC * SS(J) CALL SVP(QSAT(J),ESAT,P,RN(J)) RH(J) = QN(J) / QSAT(J) IF(RH(J).LE.0.1.OR.J.EQ.1)THEN ZSP(J)=-1. PLCL(J)=0. ELSE CALL LCL(RN(J),P,QN(J),PLCL(J),TCL(J)) CALL SP(SS,PSFC,PLCL(J),ZZ,J,MZ,ZSP(J),nlevd) ENDIF CLDLQD(J) = 0. PTNLQD(J) = 0. 200 CONTINUE DO 210 J = 2, JPBL ZLCL = ZSP(J) IF(ZSP(J).GT.0.) GOTO 220 210 CONTINUE 220 CONTINUE P2 = PSFC * SS(2) IJ = 0 IF(ZSP(2).GT.0.) THEN CALL SVP(QLCL,ELCL,PLCL(2),TCL(2)) IJ = 1 ENDIF RMXLQD = 0. DO 230 J = 2, JPBL IF(ZSP(J).GT.0.) THEN IF(PLCL(J).GT.PPBL) THEN CALL SVP(QCL,ECL,PLCL(J),TCL(J)) CALL TOPQ(QTOP,TTOP,PPBL,QCL,TCL(J),PLCL(J)) CLDLQD(J)=.5 * (QCL - QTOP) ELSE CLDLQD(J) = 0. ENDIF IF(ZZ(J).GT.ZLCL) THEN CALL SVP(QCL,ECL,PSFC*SS(J),RN(J)) CALL TOPQ(QTOP,TTOP,PPBL,QCL,RN(J),PSFC*SS(J)) PTNLQD(J) = .5 * (QCL - QTOP) * (HPBL - ZZ(J)) CLDLQD(J) = CLDLQD(J) * (HPBL - ZZ(J)) ELSE CLDLQD(J) = CLDLQD(J) * (HPBL - ZLCL) RMXLQD = AMAX1(RMXLQD,CLDLQD(J)) ENDIF ENDIF 230 CONTINUE SUMZ = 0. SUMLQD = 0. SUMPQD = 0. DO 240 J = 2, JPBL IF(J.EQ.JPBL) THEN DZCLD = HPBL - ZZ(J) + .5 * (ZZ(J) - ZZ(J-1)) ELSE DZCLD = .5 * (ZZ(J+1) - ZZ(J-1)) ENDIF SUMZ = SUMZ + DZCLD IF(ZSP(J).GT.0.) THEN SUMLQD = SUMLQD + CLDLQD(J) * DZCLD ENDIF IF(ZZ(J).GT.ZLCL.AND.ZSP(J).GT.0.) THEN SUMPQD = SUMPQD + PTNLQD(J) * DZCLD ELSE SUMPQD = SUMPQD + RMXLQD * DZCLD PTNLQD(J) = RMXLQD ENDIF 240 CONTINUE IF(SUMZ.GT.0.) THEN AVGLQD = SUMLQD / SUMZ AVGPQD = SUMPQD / SUMZ ELSE AVGLQD = 0. AVGPQD = -1. ENDIF DENS = PSFC / (287. * T2 * (1. + .61 * W2)) DO 310 J=MZBL,1,-1 P=PSFC*SS(J) PHFL=HFL(J)*1004.5*DENS PUFL=UFL(J) PVFL=VFL(J) PWFL(J)=WFL(J)*2.5E6*DENS CLWFL=0.0 ENWFL=PWFL(J) IF((PBLK(J).NE.0.0).AND.(J.NE.1)) THEN DQDZ=PWFL(J)/PBLK(J)/PR CLWFL=CLOUDK(J)*DQDZ*PR ENWFL=(PBLK(J)-CLOUDK(J))*DQDZ*PR ENDIF Q1 = TFLX(J) * 86400. Q2 = -QFLX(J) * 2.5E6 * 86400. / 1004.5 TC=RN(J) TH = RN(J)*SSK(J) QQ = QN(J) C IF(J.EQ.1) THEN C TC = RN(J+1) * SS(J+1) **(-0.286) C QQ = QN(J+1) C ENDIF CALL TDEW(QQ,P,TC,TD) IF (ZZ(J).GT.1.0E10) WRITE(*,*)' J,ZZ(J)=',J,ZZ(J) WRITE (LUNB) ZZ(J),UN(J),VN(J),TH-273.16,TC-273.16,QQ*1000., . P,SS(J),RH(J)*100.,PUFL*100.,PVFL*100.,PHFL,PWFL(J), . Q1,Q2,TD-273.16,ZSP(J),PBLK(J) 310 CONTINUE USTAR=SQRT(SQRT(TAOX*TAOX+TAOY*TAOY)) RHO=PSFC/(287.*T2*(1.+0.61*W2)) H=HEAT*1004.5*RHO E=2.5E6*EVAP*RHO HV=(HEAT+0.61*TS*EVAP)*1004.5*RHO EP1=2.5E6*EP FUP=5.67E-8*TS**4 SNOFLX = FLX1+FLX2+FLX3 SUM=H+SSOIL+E-FDOWN+FUP+SNOFLX FNET = FDOWN-FUP QM = W2 THM = TH2 QBAR = QM*ZZ2 THBAR = THM*ZZ2 IF (JPBL.GE.3) THEN DO 315 J = 3,JPBL DZL = ZZ(J)-ZZ(J-1) THVAL = RN(J)*SSK(J) QVAL = QN(J) THBAR = THBAR + 0.5*(THM+THVAL)*DZL QBAR = QBAR + 0.5*(QM+QVAL)*DZL THM = THVAL QM = QVAL 315 CONTINUE ENDIF THBAR = THBAR/ZZ(JPBL) QBAR = QBAR/ZZ(JPBL) WRITE (LUNB) FDOWN,H,SSOIL,E,FUP,SNOFLX,SUM,FNET,EP1,THBAR, . QBAR*1000.,ETOT*100. WSTAR=SIGN(1.,HEAT)*(9.806/TS*HPBL*AMAX1(1.E-15,ABS(HEAT)))**ONE3D WSC=USTAR * (1. + .47 * HPBL / XL) IF(XL.LT.0.)WSC=USTAR*(1.-0.7*HPBL/XL)**ONE3D GEOS=UG(1)*UG(1)+VG(1)*VG(1) WW2=UN(2)*UN(2)+VN(2)*VN(2) SINA=(-UN(2)*VG(1)+VN(2)*UG(1))/SQRT(GEOS*WW2) GEOS=SQRT(GEOS) XMY=0.4*USTAR/(XL*ABS(FH)) IF(Z0.NE.0.0) RO=GEOS/(ABS(FH)*Z0) SINA=FH/ABS(FH)*SINA A=0. B=0. IF(USTAR.EQ.0.)GOTO 400 A=0.4*GEOS/USTAR*SINA IF(Z0.NE.0.0) THEN B=ALOG(USTAR/GEOS*RO)-0.4*GEOS/USTAR*SQRT(1.-SINA*SINA) RO=ALOG(RO) ENDIF 400 CONTINUE TSFC = RN(1) - 273.16 HEATV = HEAT + 0.61*TS*EVAP THSTAR = -HEAT/USTAR C ---------------------------------------------------------------------- C AIR TEMPERATURE AT 2 METERS IS ESTIMATED USING THE METHOD OF C HOLTSLAG, 1988. C THE WIND SPEED AT 2 METERS IS ESTIMATED BY USING A LOGARITHMIC C WIND PROFILE WITH STABILITY CORRECTION (ALSO BY HOLTSLAG). C ---------------------------------------------------------------------- TAIR = T2M(ZZ2,TH2,TSV,THSTAR,USTAR,HEATV) - 273.16 VK = 0.4 Z0C = AMAX1(Z0,0.001) ANEM = USTAR/VK * (ALOG(2./Z0C) - FPSIM(2./XL) + FPSIM(Z01/XL)) WRITE (LUNB) TAOX,TAOY,USTAR,WSTAR,UG(1),VG(1),GEOS,SINA, . ANEM,WSC,CLC WRITE (LUNB) HPBL,XL,RIB,RIF,RIR,TSFC,TAIR,CM,CH,CG,THERM WRITE (LUNB) PRCP*3.6E6,DEW*3.6E6,ACCP*1000.,ACCD*1000. WRITE (LUNB) CMC*1000.,ESD,PC,RC,SOLARN DO 500 J=1,NSOIL TSS = TSOIL(J) - 273.16 WRITE (LUNB) DSOIL(J),WSOIL(J),TSS 500 CONTINUE WRITE (*,*) 'H,LE,Rn,G,SUM',H,E,FNET,SSOIL,SUM WRITE (*,*) 'USTAR,CLC,HBPL=',USTAR,CLC,HPBL THSFC = TSFC*SSK(1) QSFC = QN(1)*1000. T20 = RN(2)-273.16 TH20 = RN(2)*SSK(2)-273.16 Q20 = QN(2)*1000. WS20 = (UN(2)**2+VN(2)**2)**0.5 T100 = RN(6)-273.16 TH100 = RN(6)*SSK(6)-273.16 Q100 = QN(6)*1000. WS100 = (UN(6)**2+VN(6)**2)**0.5 WRITE (*,*) . '........................................................' WRITE (*,*) ' 0m T,TH,q =',TSFC,THSFC,QSFC WRITE (*,*) ' 20m T,TH,q,WS=',T20,TH20,Q20,WS20 WRITE (*,*) '100m T,Th,q,WS=',T100,TH100,Q100,WS100 WRITE (*,*) . '<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<' WRITE (*,*) return end c -------------------------------------------------------------------- * SUBROUTINE PBLL * * Model for turbulent mixing in the PBL and free atmosphere. Input * variables carried by variables in the commom block /PBL/ is mean * variables at time level "-". On return, those variables carry * tendency terms (dU/dt, etc) due to turbulent mixing. c -------------------------------------------------------------------- * MODIFICATIONS: c * 11-MAR-88 P. RUSCHER * INCORPORATE KONDO STYLE EDDY DIFFUSIVITY FOR VERY STABLE LAYER * IFCLD = CLOUD DIFFUSIVITY FLAG 0 : NO CLOUD DIFFUSION (THIS * VERSION), 1 : CLOUD DIFFUSION (FUTURE VERSION). * NOTE: CLOUD MODEL IS ACTIVE FOR EITHER FLAG; ONLY * CLOUD DIFFUSIVIRY IS TURNED OFF FOR IFCLD=0. c * 11-MAR-88 P. RUSCHER * CALCULATE A FLUX RICHARDSON NUMBER FOR USE IN SPECIFICATION OF * BOUNDARY LAYER HEIGHT. c * 2 March 1990 Jinwon Kim for the use in the version OSU2_0.FOR * * Structural change (should not alter the result, but slight * difference in the results has been caused by changing * the order of multiplication and exponentiation. Difference * is negligible.) * * Use the subroutine ROSR12 to invert tri-diagonal matrix. * * All the net fluxes above the PBL (in the free atmosphere) * will be set to be zero (to eliminate noise due to matrix * inversion.) * * This version is still PBL-only. c * 7 March 1990 Jinwon Kim * For the version: OSU2_0_2.FOR * * Zero setting the flux convergence above the PBL is now * eliminated. This results in noises in the converted * flux values above the PBL, but the magnitude is large. * * The noise above PBL may be introduced from the noise in * matrix inverting process, and revealed by vertically * integrating flux convergence to get flux at each level * in the subroutine PRINT. c * 7 March 1990 Jinwon Kim * For the version: OSU2_1.FOR * * Diffusion in the free atmosphere is added. To calculate * eddy diffusivities, one-liner formulation is used with * temporary coefficients. This subroutine will invert * matrix once for (free atmosphere+PBL) calculation. * * In the free atmosphere, PBLK(i) = Kh, and ENVK(i) = Km * * When free atmospheric diffusion is not calculated, set * any nonzero flux divergence above the PBL top to be zero. c * 27 June 1990 Jinwon Kim * * Turbulent length scales and the Prandtl number in the * stable free atmosphere derived from the CABLE 5 and 24 * May 1983 and SESAME 6 May 1979 observations. c * 1 Aug 1990 Jinwon Kim * * In the free atmosphere, when Ri<0: Use JFL formulation for * unstable case. c * * Asymptotic mixing length = 52.5 meters c * * Nonzero K at the PBL top c -------------------------------------------------------------------- * COMMON BLOCKS/ELEMENTS USED: c ------------------------------------- * /PBL/ ZZ, WINIT, MZBL, DTBL, RICR, PINK, DTNW, KOOL, IFCLD * /SFCL/ TS, WS, CM, CH, U2, V2, T2, TH2, W2, ZZ2, CBAGK * /AUXIL/ PSFC * /RAD/ CLC * /GRID/ SS, SSK * /AUX12/ PR c * OUTPUT: * /PBL/ UN, VN, RN, WN, HPBL, PBLK, CLOUDK, ZLCL, CLTOP, DCZ * /SFCL/ RIF * /AUXIL/ HEAT, TAOX, TAOY, EVAP, XL, CGH, THERM, UFLX, VFLX, * TFLX, QFLX * /HCALC/ MZH c ---------------------------------------------------------------------- * SUBROUTINES CALLED: CLOUD c ---------------------------------------------------------------------- subroutine pbll(ifree,envk,ifkpbl,rkfac,rb2) c ---------------------------------------------------------------------- parameter(nlevd=100) common /pbl/ UN(NLEVD),VN(NLEVD),RN(NLEVD),WN(NLEVD),VERT(NLEVD), . ZZ(NLEVD),WINIT(NLEVD),HPBL,MZBL,DTBL,RICR,PINK,DTNW, . KOOL,PBLK(NLEVD),CLOUDK(NLEVD),ZLCL,CLTOP,DCZ, . IFCLD common /sfcl/ TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,TSNOW,Z01,CBAGK common /auxil/ PSFC,HEAT,TAOX,TAOY,EVAP,XL,CGH,THERM,DELTAT,FH, . PRCP,DEW,ACCP,ACCD,UFLX(NLEVD),VFLX(NLEVD), . TFLX(NLEVD),QFLX(NLEVD) common /rad/ FDOWN,MO,DY,TSUN,CLC,TREF,SLO,SLA,ALBEDO,IFCRI,IFPR, . SOLARN common /grid/ ss(nlevd),ssk(nlevd),mz,mzm,noofgr common /auxi2/ ep,etot,pr common /xlun/ lunc,lunb,lund common /hcalc/ mzh,timeis real qu(nlevd),qu1(nlevd),qt(nlevd),uu(nlevd),uv(nlevd) real ut(nlevd),uw(nlevd),cool(nlevd),dry(nlevd),envk(nlevd) real dz(nlevd),dzsq(nlevd) real au(nlevd),ad(nlevd),al(nlevd) real bu(nlevd),bd(nlevd),bl(nlevd) real cu(nlevd),cd(nlevd),cl(nlevd) real yu(nlevd),yv(nlevd),work(nlevd),uut(nlevd),uvt(nlevd) real kmax c --------------------------------------------------------------- data vk,g,fak/0.4,9.806,8.5/ data betam,betah,betaz,sffrac/15.0,15.0,4.7,0.1/ data cccrt/40.0/ data c4,c5/123.6,34./ data fbc/0.5/ data c1hi,c2h,c3h,c4h/50.,-8.5,0.15,3./ data c1h0i/52.5/ data cpr0,cprs,cprus/1.5,3.08,0./ c ----------------------------------------------------------- c6=(cccrt-c5)/(100.-c5)/cccrt onet=1./3. c1h=c1hi*rkfac c1h0=c1h0i*rkfac c1hsq=c1h**2 c1h0sq=c1h0**2 c2h2=2.*c2h c ------------------------------------------------------------ * "mzblm" is the level immediately below the PBL top (HPBL) * c ------------------------------------------------------------ mzblm=mzbl-1 dcz=0. zlcl=zz(2) cltop=zz(2) bgk=1. bgk=bgk*cbagk ccon=fak*sffrac*vk binm=betam*sffrac binh=betah*sffrac cgh=0. cgw=0. therm=0. xbc=fbc*dtbl ybc=dtbl-xbc im=1 do10 i=2,mz dz(im)=zz(i)-zz(im) dzsq(im)=dz(im)*dz(im) im=i 10 continue c ---------------------------------------------------- * Compute Newtonian cooling rate based on TREF and * * DTNW. (by default, KOOL=0. in this version) * c ---------------------------------------------------- do100 j=2,mzbl qu(j)=0. qu1(j)=0. qt(j)=0. uu(j)=0. uv(j)=0. ut(j)=0. uw(j)=0. c pblk(j)=0. envk(j)=0. cloudk(j)=0. cool(j)=(rn(j)/ssk(j)-tref)*float(kool)/(10.*dtnw) dry(j)=(wn(j)-winit(j))*float(kool)/dtnw 100 continue c ----------------------------------------------------------- * Compute surface fluxes and virtual surface temperature. * c ----------------------------------------------------------- THETAS = TS * SSK(1) c dth2=th2-ts dth2=th2-thetas dw2=w2-ws taox=-cm*u2 taoy=-cm*v2 heat=-ch*dth2 evap=-ch*dw2 c heatv=heat+0.61*ts*evap heatv=heat+0.61*thetas*evap c tsv=ts*(1.+0.61*w2) c thsv=thetas*(1.+0.61*ws) thsv=thetas*(1.+0.61*w2) thl=rn(2) ustar2=sqrt(taox**2+taoy**2)+1.e-10 ustar=sqrt(ustar2) c --------------------------------- * Compute Monin-Obukhov length. * c --------------------------------- if(heatv.le.0.)then c xl=-tsv*ustar2*ustar/(g*vk*(heatv-1.e-10)) xl=-thsv*ustar2*ustar/(g*vk*(heatv-1.e-10)) else c xl=-tsv*ustar2*ustar/(g*vk*(heatv+1.e-10)) xl=-thsv*ustar2*ustar/(g*vk*(heatv+1.e-10)) endif c --------------------------------------------------- * Start calculation of boundary layer height. * c * Determine first model level virtual temperature * * and wind speed, and flux Richardson number. * c --------------------------------------------------- c t2v=t2*(1.+0.61*w2) th2v=th2*(1.+0.61*w2) vv2=u2*u2+v2*v2 vv=sqrt(vv2) c ccr=g*zz2/(t2*vv) c ccr=g*zz2/(th2*vv) ccr=g*zz2/(th2v*vv) rb2=ccr*dth2 rif=-ccr*heatv*pr/ustar2 c rif=-g*heatv*pr*zz2/(t2*ustar2*vv) ril=rif c --------------------------------------------------- * Use flux Richardson number for the lowest layer. * c --------------------------------------------------- if(rif.gt.ricr)then hpbl=zz2 jpbl=2 else jm=2 j=3 do200 jj=j,mzbl vvj=un(jj)*un(jj)+vn(jj)*vn(jj) tjv=rn(jj)*(1.+0.61*wn(jj)) c rih=g*(tjv-t2v)*zz(jj)/(thl*vvj) rih=g*(tjv-th2v)*zz(jj)/(thl*vvj) if(rih.ge.ricr)then hpbl=zz(jm)+(zz(jj)-zz(jm))*(ricr-ril)/(rih-ril) jpbl=jj goto 225 endif ril=rih jm=jj 200 continue hpbl=zz(mzbl) jpbl=mzbl endif 225 continue c --------------------------------------------------------------- * Determine more boundary layer characteristics depending on * * the stability. Stability is defined in terms of virtual * * heat flux "heatv" from the ground surface. PBL is taken to * * be stable if "heatv" <= 0,and unstable if "heatv" > 0. * c --------------------------------------------------------------- if(heatv.le.0.)then c ............................................................. * Determine fractional cloud cover in stable boundary layer. * c * If cloud flag (IFCRI) = 0 or -1, skip cloud calculation. * c ............................................................. if(ifcri.eq.1) call cloud(mzbl,clc) prinv=1./pr fak1=ustar*hpbl*vk*dtbl else c ................................................... * Unstable case: Compute velocity scale, thermal * * temperature excess and counter-gradient terms. * c ................................................... fmt=(1.-binm*hpbl/xl)**onet fht=sqrt(1.-binh*hpbl/xl)/pr wsc=ustar*fmt therm=fak/(hpbl*wsc) cgh=heat*therm cgw=evap*therm c ................................................... * Put cap on the counter-gradient factor for heat * * (Deardorff, 1966) and turn off CGW. * c cgw=0.0 *after* therm=... in the section below?; c cgw important for the 'therm' term below? c --- c 12-May-1993, again trun on the 'countergradient' (nonlocal) mixing c of moisture, Holtslag and Ek, 1993. c ................................................... cgh=amin1(cgh,6.5e-4) c cgw=0. c therm=(cgh+0.61*ts*cgw)*hpbl therm=(cgh+0.61*thetas*cgw)*hpbl c cgw=0. c ................................................................... * Improve HPBL estimation under (unstable) convective conditions. * c ................................................................... c ril=-g*therm*zz2/(tsv*vv2) ril=-g*therm*zz2/(thsv*vv2) c tlv=t2v+therm tlv=th2v+therm jm=2 j=3 do340 jj=j,mzbl vv2=un(jj)*un(jj)+vn(jj)*vn(jj) tjv=rn(jj)*(1.+0.61*wn(jj)) c rih=g*(tjv-tlv)*zz(jj)/(vv2*tsv) rih=g*(tjv-tlv)*zz(jj)/(vv2*thsv) if(rih.ge.ricr)then hpbl=zz(jm)+(ricr-ril)/(rih-ril)*dz(jm) jpbl=jj goto 345 endif ril=rih jm=jj 340 continue hpbl=zz(mzbl) jpbl=mzbl 345 continue c ............................................................... * Determine fractional cloud cover under unstable conditions. * * If cloud flag (IFCRI) = 0 or -1, skip cloud calculation * c ............................................................... if(ifcri.eq.1) call cloud(mzbl,clc) kmax=0. if(clc.ge.0.)then kmax=c4*(clc-c5)/(100.-c5)*float(ifclc) else kmax=c4*c6*clc*float(ifcld) endif c .............................................. * Compute integrals of diffusivity and store * * temporarily in QU work array. * c .............................................. vkht=hpbl*vk*dtbl fak1=ustar*vkht fak2=wsc*vkht prinv=1./(fmt/fht+ccon) endif c --------------------------------------------------------- * Define "MZH" as grid level at/immediately above HPBL. * c * Calculate K's upto MZH <= MZBL. * c --------------------------------------------------------- zp=zz2 mzh=jpbl mzhp=mzh+1 mzhm=mzh-1 if(kmax.gt.0.)then zhb=zlcl/hpbl center=(9.+5.*zhb)/14. sigsq=(1.-center)**2 endif c -------------------------------------------------------------------- * set pblk array to zero. c -------------------------------------------------------------------- do401 i=2,mzbl pblk(i)=0. 401 continue c ---------------------------------------------- * Calculate eddy diffusivities within the PBL * c ---------------------------------------------- do400 j=2,mzblm jp=j+1 zzh=0. zm=zp zp=zz(jp) du=un(jp)-un(j) dv=vn(jp)-vn(j) dt=rn(jp)-rn(j) tbr=.5*(rn(jp)+rn(j)) uvsq=du*du+dv*dv+1.e-10 ril=g*dt*dz(j)/(tbr*uvsq) uvabs=sqrt(uvsq) if((kmax.eq.0.).and.(bgk.eq.0.).and.(zp.gt.hpbl))zp=hpbl * Skip below if no clouds, no background diffusivity, * * and above the boundary layer. * if(.not.(kmax.eq.0..and.bgk.eq.0..and.zm.ge.hpbl))then z=0.5*(zm+zp) zh=z/hpbl zl=z/xl if(zh.le.1.)then zzh=1.-zh zzh=zzh**pink endif if(heatv.le.0.)then c ............................................................... * if the flag "ifkpbl"=1: * determine Km,h for stable case by using PBL parameterization * * and oneliner model. Then use larger K for calculation. * c * in a very stable situation, use formula by Ruscher (1990) * * with a Condo modification (1978), that is, for z/L > 1, * * set z/L = 1. (PBL model) * c * for oneliner model in the PBL, limit the mixing length by * * "kz" where k is the von Karman constant. * c * use max. of oneliner K or PBL model K above z/h = 0.3 * c ............................................................... pktmp=0. ektmp=0. c .... K by free atm. diffusion model ........... if(ifkpbl.eq.1)then xlh1=amin1(vk*z,c1h) if(ril.le.0.)then prdt=cpr0 ctemp=1. else prdt=cpr0+cprs*ril ctemp=(exp(c2h*ril)+c3h/(ril+c4h)) endif xlh=xlh1*ctemp xlh=amin1(xlh,c1h) xlhsq=xlh**2 pktmp=xlhsq*uvabs/dz(j) ektmp=pktmp*prdt endif c .... K by PBL parameterization .............. phim=1.+betaz*amax1(zl,1.) pblk(jp)=fak1*zh*zzh/(phim*dtbl) c .... select K above the level z=0.3*zh ...... if(zh.gt.0.3.and.ifkpbl.eq.1)then pblk(jp)=amax1(pblk(jp),pktmp) pblk(jp)=amax1(pblk(jp),bgk) endif envk(jp)=amax1(pblk(jp),ektmp) qu(jp)=pblk(jp)*dz(j)*dtbl qu1(jp)=qu(jp) else c .................................................... * Determine eddy diffusivities for unstable case. * * If clouds are present, include cloud diffusivity. * c .................................................... if(zh.lt.sffrac)then * Eddy diffusivity within the surface layer * if((kmax.gt.0.).and.(z.ge.zlcl))then cloudk(jp)=kmax*exp(-(zh-center)**2/sigsq) endif phim=1./(1.-betam*zl)**onet pblk(jp)=fak1*zh*zzh/(phim*dtbl) envk(jp)=pblk(jp) zt=dz(j)*dtbl qu(jp)=pblk(jp)*zt qu1(jp)=envk(jp)*zt else * Eddy diffusivity within the mixed layer. * if((kmax.gt.0.).and.(z.ge.zlcl).and.(z.le.cltop))then cloudk(jp)=kmax*exp(-(zh-center)**2/sigsq) endif crtk=fak2*zh*zzh/dtbl if(zh.le.1.)then if(zh.gt.0.3)crtk=amax1(crtk,bgk) else crtk=bgk endif pblk(jp)=crtk+cloudk(jp) envk(jp)=crtk zt=dz(j)*dtbl qu(jp)=pblk(jp)*zt qu1(jp)=envk(jp)*zt endif endif endif 400 continue c ------------------------------------------------------ * Calculate eddy diffusivities in the free atmosphere * * using one-liner model. To run this part, set the * * variable "IFREE" to be 1 in the data statement. * * Otherwise, it will be skipped. Leave K=0 at the * * top of PBL. * c * Note that mixing length formulation is derived only * * for stable stratification. For unstable case, the * * form of mixing length used in here has no base to * * be justified (it is arbitrary with cap on). * c ------------------------------------------------------ if(ifree.eq.1)then im=mzh do410 i=mzhp,mzblm du=un(i)-un(im) dv=vn(i)-vn(im) dt=rn(i)-rn(im) tbr=.5*(rn(i)+rn(im)) uvsq=du*du+dv*dv+1.e-10 ril=g*dt*dz(im)/(tbr*uvsq) uvabs=sqrt(uvsq) if(ril.gt.0.)then ! Stable condition prdt=cpr0+cprs*ril xlh=c1h*(exp(c2h*ril)+c3h/(ril+c4h)) xlhsq=xlh**2 pblk(i)=xlhsq*uvabs/dz(im) envk(i)=pblk(i)*prdt else ! Unstable free atm. Use JFL dzz=2.*dz(im)/(zz(i)+zz(im)) czsq=(c1h0*dz(im))**2*sqrt(dzz)*((1.+dzz)**onet-1.)**1.5 ff=ril/(1.+75.*czsq*sqrt(abs(ril))) fht=1.-15.*ff fmt=1.-10.*ff coeff=c1h0sq*uvabs/dz(im) pblk(i)=coeff*fht envk(i)=coeff*fmt c write(*,*)'Ri<0 in f.a. L=',i,pblk(i),envk(i) endif zt=dz(im)*dtbl qu(i)=pblk(i)*zt qu1(i)=envk(i)*zt im=i 410 continue endif * Set no flux condition at the top level (mzbl) * pblk(mzbl)=0. envk(mzbl)=0. qu(mzbl)=0. qu1(mzbl)=0. c -------------------------------------------------- * Subtract Ts and CGW value from pot. temp. and q * c -------------------------------------------------- if(heatv.le.0.)then do420 j=2,mzbl c rn(j)=rn(j)-ts rn(j)=rn(j)-thetas wn(j)=wn(j)-ws 420 continue else do425 j=2,mzbl dz2=zz(j)-zz2 c rn(j)=rn(j)-ts-cgh*dz2 rn(j)=rn(j)-thetas-cgh*dz2 wn(j)=wn(j)-ws-cgw*dz2 425 continue endif c-------------------------------------------------------- * Solve diffusion equation. Use the subroutine ROSR12. * * In the following, trailing letters (u,d,l) following * * coefficients (a,b,c) denote the position of vars. * * in a coefficient matrix (u=upper diagonal, d=diag.,* * and l=lower diagonsl). * c * The equation to solve is: * * CX=Y where C=A+B is the coefficient matrix. * c * Note: The following procedure is only for PBL calc. * * All the coefficients in the free atmosphere are set* * to be zero, and any noise in the free atmosphere * * as a results of matrix inversion is set to be zero * * after inverting the matrix. * c ------------------------------------------------------- * First solve for momentum (u, and v wind). * c ------------------------------------------------------- * bu(j)=b(j,j+1) * do500 j=2,mzblm jp=j+1 bu(j)=-qu1(jp)/dzsq(j) au(j)=dz(j)/6. 500 continue * Assign UBC and LBC. (no flux at the top) * bu(mzbl)=0. au(mzbl)=0. bl(2)=0. al(2)=0. * bl(j)=b(j,j-1): the matrix is symmetric * jm=2 do505 j=3,mzbl bl(j)=bu(jm) al(j)=au(jm) jm=j 505 continue * The first and last diagonal components * bd(2)=-bu(2) bd(mzbl)=-bl(mzbl) ad(2)=2.*au(2) ad(mzbl)=2.*al(mzbl) * Diagonal components of b(): bd(j)=b(j,j) * do510 j=3,mzblm bd(j)=-(bu(j)+bl(j)) ad(j)=2.*(au(j)+al(j)) 510 continue c ---------------------------------------------------------- * The right-hand-side of the elementary equation and * * the boundary conditions are specified here. The LBC * * is a flux condition (Natural, Neumann), and is not * * matched exactly. The top B.C. is a zero flux condition.* c ---------------------------------------------------------- yu(1)=ad(2)*un(2)+au(2)*un(3)+taox*dtbl yv(1)=ad(2)*vn(2)+au(2)*vn(3)+taoy*dtbl yu(mzblm)=al(mzbl)*un(mzblm)+ad(mzbl)*un(mzbl) yv(mzblm)=al(mzbl)*vn(mzblm)+ad(mzbl)*vn(mzbl) jm=2 do515 j=3,mzblm jp=j+1 yu(jm)=al(j)*un(jm)+ad(j)*un(j)+au(j)*un(jp) yv(jm)=al(j)*vn(jm)+ad(j)*vn(j)+au(j)*vn(jp) jm=j 515 continue * The left-hand-side matrix is A+B * jm=1 do520 j=2,mzbl cu(jm)=bu(j)+au(j) cd(jm)=bd(j)+ad(j) cl(jm)=bl(j)+al(j) jm=j 520 continue * Call the subroutine "ROSR12" to invert the C=A+B matrix with Y * call rosr12(uut,cl,cd,cu,yu,work,mzblm) call rosr12(uvt,cl,cd,cu,yv,work,mzblm) * Calculate dU/dt and dV/dt by diffusive process. dU/dt * * and dV/dt are set to be zero at the top and bottom levels. * un(1)=0. vn(1)=0. un(mzbl)=0. vn(mzbl)=0. do525 j=1,mzblm uu(j)=uut(j) uv(j)=uvt(j) 525 continue c ------------------------------------------------------------------- * Solve for potential temperature. Array a's are reused (note that * * al(),ad(),au() are only function of grid spacing, hence they * * are unchanged for all variables. Only the coefficients b's and * * y's are needed to be recalculated. * c * The procedures are identical for momentum flux except that eddy * * diffusivities are for heat. * c ------------------------------------------------------------------- do530 j=2,mzh jp=j+1 c bu(j)=-qu1(jp)*prinv/dzsq(j) bu(j)=-qu(jp)*prinv/dzsq(j) 530 continue do531 j=mzhp,mzblm jp=j+1 bu(j)=-qu(jp)/dzsq(j) 531 continue bu(mzbl)=0. bl(2)=0. jm=2 do535 j=3,mzbl bl(j)=bu(jm) jm=j 535 continue * The first and last diagonal elements * bd(2)=-bu(2) bd(mzbl)=-bl(mzbl) do540 j=3,mzblm bd(j)=-(bu(j)+bl(j)) 540 continue * Assign the right-hand-side of equation with B.C.'s * yu(1)=ad(2)*rn(2)+au(2)*rn(3)+heat*dtbl yu(mzblm)=al(mzbl)*rn(mzblm)+ad(mzbl)*rn(mzbl) jm=2 do545 j=3,mzblm jp=j+1 yu(jm)=al(j)*rn(jm)+ad(j)*rn(j)+au(j)*rn(jp) jm=j 545 continue * L.H.S. = A+B * jm=1 do550 j=2,mzbl cu(jm)=bu(j)+au(j) cd(jm)=bd(j)+ad(j) cl(jm)=bl(j)+al(j) jm=j 550 continue * Invert the matrix * call rosr12(uut,cl,cd,cu,yu,work,mzblm) rn(1)=0. rn(mzbl)=-cool(mzbl) do555 j=1,mzblm ut(j)=uut(j) 555 continue c -------------------------------------------------------------------- * Solve for moisture. Procedure is the same. Only K's are different * c -------------------------------------------------------------------- do560 j=2,mzh jp=j+1 bu(j)=-qu(jp)*prinv/dzsq(j) 560 continue do561 j=mzhp,mzblm jp=j+1 bu(j)=-qu(jp)/dzsq(j) 561 continue bu(mzbl)=0. * Assign "bl"'s * bl(2)=0. jm=2 do565 j=3,mzbl bl(j)=bu(jm) jm=j 565 continue * Assign "bd"'s * bd(2)=-bu(2) bd(mzbl)=-bl(mzbl) do570 j=3,mzblm bd(j)=-(bu(j)+bl(j)) 570 continue * Assign L.H.S. of the equation * jm=1 do575 j=2,mzbl cu(jm)=bu(j)+au(j) cd(jm)=bd(j)+ad(j) cl(jm)=bl(j)+al(j) jm=j 575 continue * Assign R.H.S. of the equation, together with B.C.'s * yu(1)=ad(2)*wn(2)+au(2)*wn(3)+evap*dtbl yu(mzblm)=al(mzbl)*wn(mzblm)+ad(mzbl)*wn(mzbl) jm=2 do580 j=3,mzblm jp=j+1 yu(jm)=al(j)*wn(jm)+ad(j)*wn(j)+au(j)*wn(jp) jm=j 580 continue * Invert the coefficient matrix with Y (R.H.S.) * call rosr12(uut,cl,cd,cu,yu,work,mzblm) wn(1)=0. wn(mzbl)=-dry(mzbl) do585 j=1,mzblm uw(j)=uut(j) 585 continue c ----------------------------------------------- * Calculate flux divergence term for each level * It will be arbitrarily set to be zero above * PBL (as was done in Ruscher's version.) c ----------------------------------------------- im=1 do620 i=2,mzblm un(i)=(uu(im)-un(i))/dtbl vn(i)=(uv(im)-vn(i))/dtbl rn(i)=(ut(im)-rn(i))/dtbl-cool(i) wn(i)=(uw(im)-wn(i))/dtbl-dry(i) im=i 620 continue if(ifree.ne.1)then do625 i=mzhp,mzblm un(i)=0. vn(i)=0. rn(i)=-cool(i) wn(i)=-dry(i) 625 continue endif un(mzbl)=0. vn(mzbl)=0. rn(mzbl)=-cool(mzbl) wn(mzbl)=-dry(mzbl) c ------------------------- * End of flux calculation c ------------------------- uflx(1)=0. vflx(1)=0. tflx(1)=-cool(1) qflx(1)=-dry(1) do650 j=2,mzbl uflx(j)=un(j) vflx(j)=vn(j) tflx(j)=rn(j)+cool(j) qflx(j)=wn(j)+dry(j) 650 continue c ------------------------------------- return end C ********************************************************************** C SUBROUTINE SUN C -------------- C TOTAL DOWNWARD RADIATION DETERMINED FROM NET INCOMING SOLAR C RADIATION, INCOMING SOLAR * (1 - ALBEDO), PLUS DOWNWARD C TERRESTRIAL RADIATION, BOTH MODIFIED BY CLOUDS. C ---------------------------------------------------------------------- C SOURCE: 6-FEB-87, REPLACES PREVIOUS SUN SUBROUTINE. C 2-AUG-88, P. RUSCHER C MODIFIED TO TAKE CLOUD FEEDBACK OUT OF THIS SUBROUTINE. C 18-MAR-88, P. RUSCHER C REFERENCE TEMPERATURE IS NOW MODEL LEVEL 3 FOR A18 GRID, C CALCULATED AS CLOSEST TO ~200 M FOR OTHERS OTHER CHANGES C COSMETIC. C 20-APR-88, P. RUSCHER C CALCULATION OF EFFECTIVE ATMOSPHERIC TEMPERATURE USED FOR C DOWNWARD LONGWAVE RADIATION CHANGED TO REFLECT VARIABLE GRID C STRUCTURE OF MODEL. C 15-MAY-89, M. EK C CLOUD-RADIATION FLAG (IFCRI) ADDED TO MODEL SO CLOUD FEEDBACK C COULD BE REMOVED FROM THIS SUBROUTINE WITHOUT SETTING CLOUD C COVER EQUAL TO ZERO. C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: NONE C OUTPUT: NONE C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /FIELDS/ RNM, QNM C /PBL/ ZZ C /GRID/ SS, MZ, NOOFGR C /RAD/ MO, DY, HO, CLC, SLO, SLA, R, IFCRI, IFPR C /SFCL/ TS C /AUXIL/ PSFC, HEAT, EVAP C OUTPUT: C /RAD/ FDOWN, TREF C ---------------------------------------------------------------------- C SUBROUTINES CALLED: NONE C ********************************************************************** SUBROUTINE SUN PARAMETER (NLEVD=100,NLEVI=100) INTEGER *4 IHA COMMON /FIELDS/ UN(NLEVD),UNM(NLEVI),VN(NLEVD),VNM(NLEVI), . RN(NLEVD),RNM(NLEVI),QN(NLEVD),QNM(NLEVI), . UG(NLEVD),VG(NLEVD) COMMON /PBL/ DUDT(NLEVD),DVDT(NLEVD),DTDT(NLEVD),DWDT(NLEVD), . VERT(NLEVD),ZZ(NLEVD),WINIT(NLEVD),HPBL,MZBL,DELT2, . RICR,PINK,DTNW,KOOL,PBLK(NLEVD),CLOUDK(NLEVD), . ZLCL,CLTOP,DCZ,IFCLD COMMON /GRID/ SS(NLEVD),SSK(NLEVD),MZ,MZM,NOOFGR COMMON /SFCL/ TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,TSNOW,Z01,CBAGK COMMON /AUXIL/ PSFC,HEAT,TAOX,TAOY,EVAP,XL,CG,THERM,DELTAT,FH, . PRCP,DEW,ACCP,ACCD,UFLX(NLEVD),VFLX(NLEVD), . TFLX(NLEVD),QFLX(NLEVD) COMMON /RAD/ FDOWN,MO,DY,HO,CLC,TREF,SLO,SLA,ALBEDO,IFCRI,IFPR, . SOLARN DIMENSION TSIN(360),TCOS(360) DATA PI,XA1,XA2 /3.14159,990.,-30./ DATA SIGMA,C2 /5.67E-8,60./ DATA DNF,HNF /.986,15./ DATA TRANSM,N /0.7,1./ C ---------------------------------------------------------------------- C 0. IF CLOUD-RADIATION INTERACTION FLAG (IFCRI) = 0 THEN THERE IS C NO CLOUD-RADIATION INTERACTION. SET FLAG = 0. C OTHERWISE SET FLAG = 1. C ---------------------------------------------------------------------- IF (IFCRI .EQ. 0) THEN FLAG = 0. ELSE FLAG = 1. ENDIF C ---------------------------------------------------------------------- C 1. COMPUTE TABLES AND CONSTANTS. C ---------------------------------------------------------------------- 100 DO 110 I=1,360 ANG = FLOAT(I-1) * PI/180. TSIN(I) = SIN(ANG) TCOS(I) = COS(ANG) 110 CONTINUE SQT3 = SQRT(3.) C ---------------------------------------------------------------------- C 2. COMPUTE DAY NUMBER. C ---------------------------------------------------------------------- 200 DN = 30.*(MO-1) + DY C ---------------------------------------------------------------------- C 3. COMPUTE SOLAR ALTITUDE, SFI = SIN(SOLAR ALTITUDE). C ---------------------------------------------------------------------- IL = MOD(IFIX((DN+9.)*DNF)+360,360)+1 DECL = -23.5*TCOS(IL) + 360. IHA = -SLO+HO*HNF + 540. XCH = TCOS(MOD(IHA,360)+1) ISLA = MOD(IFIX(SLA)+360,360) + 1 IDEC = MOD(IFIX(DECL),360) + 1 SFI = TSIN(IDEC)*TSIN(ISLA) + TCOS(IDEC)*TCOS(ISLA)*XCH C ---------------------------------------------------------------------- C 4. DETERMINE SOLAR RADIATION. C COMPUTE INCOMING SOLAR RADIATION FOR CLEAR SKY. C REDUCE INCOMING SOLAR RADIATION FOR CLEAR SKY DUE TO CLOUDS. C LET SOME SOLAR RADIATION PENETRATE THRU THE CLOUDS (=TRANSM). C SOLARN IS THE NET SOLAR RADIATION AT THE SURFACE (DOWNWARD SOLAR C LESS REFLECTED). C ---------------------------------------------------------------------- TRANSM = 0.06 + 0.17*SIN(SFI) OPAQUE = (1.0 - TRANSM) CLEAR = AMAX1(XA1*SFI+XA2,0.) c SOLAR = CLEAR * (1. - TRANSM*(CLC**N)*FLAG) SOLAR = CLEAR * (1. - OPAQUE*CLC*FLAG) SOLARN = (1. - ALBEDO) * SOLAR C ---------------------------------------------------------------------- C 5.1. DETERMINE DOWNWARD LONGWAVE RADIATION. C IF IFPR = 0 USE SPECIFIED EFFECTIVE ATMOSPHERIC RADIATING C TEMPERATURE (TREF) TO DETERMINE DOWNWARD LONGWAVE RADIATION. C ---------------------------------------------------------------------- IF (IFPR .EQ. 0) THEN TDOWN = SIGMA * TREF**4. C ---------------------------------------------------------------------- C 5.2. IF IFPR = 1 THEN USE MORE SOPHISTICATED PARAMETERIZED C RADIATION FOR LONGWAVE FOLLOWING SATTERLUND (1979) WITH CLOUD C CORRECTION BY PALTRIDGE AND PLATT (1976). C FIND MODEL LEVEL CLOSEST TO 200 METERS; THIRD MODEL LEVEL C (SECOND ABOVE SFC) FOR GRID #6. C ---------------------------------------------------------------------- ELSE IF (NOOFGR.EQ.6) THEN K2=3 ELSE KJ = 1 220 KJ = KJ+1 IF (ZZ(KJ).LT.200.) GO TO 220 IF (ZZ(KJ)-200. .GT. 200.-ZZ(KJ-1)) THEN K2 = KJ-1 ELSE K2 = KJ ENDIF ENDIF C ---------------------------------------------------------------------- C 5.3. DEFINE REFERENCE TEMPERATURE, HUMIDITY, AND PRESSURE AT C THIS LEVEL. VAPOR PRESSURE (EREF) IN MB; RESTRICT EREF ABOVE C 0.05 MB FOR COMPUTATION STABILITY. C ---------------------------------------------------------------------- XTREF = RNM(K2) QREF = QNM(K2) PRESS = PSFC * SS(K2) EREF = PRESS * QREF * 0.01/0.622 EREF = AMAX1(EREF,.05) XT2 = XTREF * XTREF T4 = XT2 * XT2 C ---------------------------------------------------------------------- C 5.4. DETERMINE EMMISSIVITY (EMSVTY) OF THE ATMOSPHERE. C MODIFY DOWNWARD LONGWAVE RADIATION DUE TO CLOUDS. C ---------------------------------------------------------------------- EMSVTY = 1.08 * (1. - EXP(-EREF ** (XTREF / 2016.))) TDOWN = EMSVTY * SIGMA * T4 + C2 * CLC * FLAG ENDIF C ---------------------------------------------------------------------- C 6. COMPUTE TOTAL DOWNWARD RADIATIVE FLUX = C TDOWN (DOWNWARD LONGWAVE) + INCOMING NET SOLAR (SOLARN). C ---------------------------------------------------------------------- FDOWN = TDOWN + SOLARN RETURN END C ********************************************************************** C SUBROUTINE CLOUD C ----------------- C THIS ROUTINE CALCULATES THE FRACTIONAL CLOUD COVER AS A FUNCTION C OF THE MAXIMUM RELATIVE HUMIDITY IN THE BOUNDARY LAYER AND THE C RELATIVE HUMIDITY VARIANCE AT THE SAME LEVEL; THE VARIANCE C DEPENDS ON BOTH TURBULENT AND MESOSCALE (SUBGRID) EFFECTS. C ---------------------------------------------------------------------- C SOURCE: EK AND MAHRT (1991) C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: JPBL C OUTPUT: CLC C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: C /EK/ JUST FOR TESTING. CAN BE DELETED HERE AND IN SUB PRINT C /SFCL/ TS, TH2, W2 C /PBL/ RN, WN, ZZ, HPBL, PBLK C /AUXIL/ PSFC, HEAT, TAOX, TAOY, EVAP C /GRID/ SS, SSK C ---------------------------------------------------------------------- C SUBROUTINES CALLED: SVP C ********************************************************************** SUBROUTINE CLOUD(JPBL,CLC) PARAMETER(NLEVD=100) COMMON /SFCL/ TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,TSNOW,Z01,CBAGK COMMON /PBL/ UN(NLEVD),VN(NLEVD),RN(NLEVD),WN(NLEVD),VERT(NLEVD), . ZZ(NLEVD),WINIT(NLEVD),HPBL,MZBL,DTBL,RICR,PINK,DTNW, . KOOL,PBLK(NLEVD),CLOUDK(NLEVD),ZLCL,CLTOP,DCZ, . IFCLD COMMON /AUXIL/ PSFC,HEAT,TAOX,TAOY,EVAP,XL,CGH,THERM,DELTAT,FH, . PRCP,DEW,ACCP,ACCD,UFLX(NLEVD),VFLX(NLEVD), . TFLX(NLEVD),QFLX(NLEVD) COMMON /GRID/ SS(NLEVD),SSK(NLEVD),MZ,MZM,NOOFGR COMMON /EK/ RHMAXX,LVLMAX C ---------------------------------------------------------------------- C AX,BX ARE REGRESSION COEFFICIENTS TO GET TURBULENT STD DEV RH FROM C MOISTURE FLUX, STD DEV W, AND QSAT. C SRHMS2 IS MESOSCALE VARIANCE OF RH. C ---------------------------------------------------------------------- DATA G,AX,BX,SRHMS2 /9.806,9.81,0.00024,0.0016/ DATA A0,A1,A3,A5,A7,A9 . /0.5,0.3938824,-5.843166E-2,6.139052E-3,-3.394664E-4,7.433019E-6/ RHMAX = 0.0 WSTAR2 = 0.0 SIGWWS = 0.0 SIGWSN = 0.0 SIGMAW = 0.0 SRHTR2 = 0.0 IFLAG = 0 MZHX = JPBL HEATV = HEAT + 0.61*TS*EVAP TH2V = TH2 * (1+0.61*W2) USTAR2 = SQRT(TAOX*TAOX + TAOY*TAOY) + 1.0E-10 C ---------------------------------------------------------------------- C DETERMINE MAXIMUM RH AND LEVEL IN BOUNDARY LAYER. FIRST TIME STEP C MZHX=0 SO RESET AT 2; AFTER THE FIRST TIME STEP MZHX IS >=2. C ---------------------------------------------------------------------- MZHX = MAX0(2,MZHX) DO 100 J=2,MZHX P = PSFC*SS(J) T = RN(J)/SSK(J) CALL SVP(QS,ES,P,T) RH = WN(J)/QS IF (RH .GT. RHMAX) THEN QSMAX = QS RHMAX = RH IMAX = J ENDIF 100 CONTINUE RHMAXX = RHMAX LVLMAX = IMAX C ---------------------------------------------------------------------- C DETERMINE MOISTURE FLUX AT IMAX (LEVEL OF RH MAXIMUM) IN ORDER TO C EVENTUALLY DETERMINE STD DEV RH AT THAT LEVEL. C ---------------------------------------------------------------------- QFLX1 = PBLK(IMAX-1) * (WN(IMAX)-WN(IMAX-1))/(ZZ(IMAX)-ZZ(IMAX-1)) QFLX2 = PBLK(IMAX) * (WN(IMAX+1)-WN(IMAX)) / (ZZ(IMAX+1)-ZZ(IMAX)) DZ1 = 0.5 * (ZZ(IMAX) + ZZ(IMAX-1)) DZ2 = 0.5 * (ZZ(IMAX+1) + ZZ(IMAX)) QFLXMX = (ZZ(IMAX)-DZ1)/(DZ2-DZ1) * QFLX2 + . (DZ2-ZZ(IMAX))/(DZ2-DZ1) * QFLX1 C ---------------------------------------------------------------------- C DETERMINE VERTICAL VELOCITY VARIANCE (SIGMAW). C SIGMAW (SIGWWS) FOR UNSTABLE CONDITIONS. C ---------------------------------------------------------------------- IF (HEATV .GT. 0.0) THEN WSTAR2 = (G*HPBL*HEATV/TH2V)**0.666 XPBLWS = 1.8*((ZZ(IMAX)/zz(mzhx))**0.666)* . (1.0-0.8*ZZ(IMAX)/zz(mzhx))**2 SIGWWS = (WSTAR2 * XPBLWS)**0.5 ENDIF C ---------------------------------------------------------------------- C SIGMAW (SIGWSN) FOR STABLE AND NEUTRAL CONDITIONS. C SIGMAW = MAXIMUM OF SIGWWS AND SIGWSN FOR UNSTABLE C CONDITIONS, OR JUST SIGWSN FOR STABLE CONDITIONS. C ---------------------------------------------------------------------- IF (ZZ(IMAX) .NE. ZZ(MZHX)) THEN XPBLSN = 2.5 * (1 - (ZZ(IMAX)/zz(mzhx))**0.6) SIGWSN = (USTAR2 * XPBLSN)**0.5 SIGMAW = AMAX1(SIGWWS,SIGWSN) XX = QFLXMX/(SIGMAW*QSMAX) C ---------------------------------------------------------------------- C STABLE OR NEUTRAL AND Z=H, THEN SIGMAW=0, SO SRHTR2=0. C ---------------------------------------------------------------------- ELSE IF (HEATV .LE. 0.0) THEN XX = 0.0 C ---------------------------------------------------------------------- C UNSTABLE AND Z=H, SIGMAW=SIGWWS. C ---------------------------------------------------------------------- ELSE XX = QFLXMX/(SIGWWS*QSMAX) ENDIF C ---------------------------------------------------------------------- C DETERMINE RH STD DEV (SRH) AT LEVEL OF RH MAXIMUM (IMAX). C TURBULENT CONTRIBUTION TO SRH DEPENDS ON MOISTURE FLUX (QFLXMX), C VERTICAL VELOCITY VARIANCE (SIGMAW), AND SATURATION SPECIFIC C HUMIDITY (QSMAX), PLUS A FIXED VALUE (BX). C MESOSCALE CONTRIBUTION DEPENDS ON HORIZONTAL SCALE (0.04**2 USED C FOR 100 KM HORIZONTAL SCALE). C ---------------------------------------------------------------------- SRHTR2 = AX*(XX**2) + BX SRH = (SRHTR2 + SRHMS2)**0.5 C ---------------------------------------------------------------------- C DETERMINE STANDARDIZED (STATISTICAL) Z VALUE FROM MEAN RH (RHMAX) C AND STD DEV OF RH (SRH). C IF Z < -3 SET CLC=0; IF Z > 3 SET CLC=1; THEN TURN ON FLAG TO SKIP C THE NUMERICAL CALCULATION OF CLC. C ---------------------------------------------------------------------- Z = (RHMAX-1.) / SRH IF (Z .LT. -3.0) THEN CLC = 0.0 IFLAG = 1 ELSE IF (Z .GT. 3.0) THEN CLC = 1.0 IFLAG = 1 ENDIF C ---------------------------------------------------------------------- C POLYNOMIAL FIT TO GAUSSIAN DISTRIBUTION. FRACTIONAL CLOUD COVER C (CLC) IS THE AREA UNDER THE NORMAL CURVE GREATER THAN Z. C ---------------------------------------------------------------------- IF (IFLAG .EQ. 0) THEN CLC = A0 + A1*Z + A3*Z**3 + A5*Z**5 + A7*Z**7 + A9*Z**9 CLC = AMAX1(0.0,CLC) CLC = AMIN1(CLC,1.0) ENDIF RETURN END C ********************************************************************** C SUBROUTINE CANRES C ----------------- C THIS ROUTINE CALCULATES THE CANOPY RESISTANCE WHICH DEPENDS C ON INCOMING SOLAR RADIATION, AIR TEMPERATURE AND ATMOSPHERIC C WATER VAPOR PRESSURE DEFICIT AT THE LOWEST MODEL LEVEL. C ---------------------------------------------------------------------- C SOURCE: JARVIS (1976), Jacquemin and Noilhan (1990 BLM) C ---------------------------------------------------------------------- C ARGUMENTS: C INPUT: RR, DELTA C ---------------------------------------------------------------------- C COMMON BLOCKS/ELEMENTS USED: C INPUT: /RAD/ SOLARN,ALBEDO C /SFCL/ CH,T2,W2 C /AUXIL/ PSFC C /GRID/ SS C OUTPUT: /SOIL1/ PC,RC C ---------------------------------------------------------------------- C SUBROUTINES CALLED: SVP C ********************************************************************** SUBROUTINE CANRES(RR,DELTA,ircflg) PARAMETER(NLEVD=100,NSOLD=10) COMMON /RAD/ FDOWN,MO,DY,HO,CLC,TREF,SLO,SLA,ALBEDO,IFCRI,IFPR, . SOLARN COMMON /SOIL/ WSOIL(NSOLD),CMC,NSOIL,SSOIL,TSOIL(NSOLD),ESD COMMON /SOIL1/ IFTC,TSO0,SIGMAF,TSOSAT,TWILT,SCANOP, . PC,RC,KWILT,CFACTR,TSOREF COMMON /SOIL2/ ZSOIL(NSOLD),ISOIL,IFSOIL COMMON /SFCL/ TS,WS,CM,CH,U2,V2,T2,TH2,W2,ZZ2,Z0,Z0H,BETA,RIR, . RIB,RIF,TSNOW,Z01,CBAGK COMMON /GRID/ SS(NLEVD),SSK(NLEVD),MZ,MZM,NOOFGR COMMON /AUXIL/ PSFC,HEAT,TAOX,TAOY,EVAP,XL,CG,THERM,DELTAT,FH, . PRCP,DEW,ACCP,ACCD,UFLX(NLEVD),VFLX(NLEVD), . TFLX(NLEVD),QFLX(NLEVD) COMMON /HCALC/ MZH,TIMEIS COMMON /PLANT/ RCMIN DIMENSION PART(10) character tab*1 C ---------------------------------------------------------------------- C FOR PINE FOREST... C FOR SOLAR FUNCTION (RCS) USE THESE VALUES BY PINTY ET AL (89 MWR). C (A1,A2,XLAI,RSMIN,RXMAX) C TLOW,TOPT,THI FROM ??? C C1 value from Noilhan and Planton (89 MWR) C ---------------------------------------------------------------------- C DATA A1,A2,XLAI,RSMIN,RSMAX C . /1.,0.0055,3.5,350.,2000./ C DATA TLOW,TOPT,THI /268.15,283.15,307.15/ C DATA C1 /0.00025/ C ---------------------------------------------------------------------- C FOR PINE FOREST USE THESE VALUES FROM JACQUEMIN AND NOILHAN (90 BLM) C ---------------------------------------------------------------------- DATA A1,A3,A4,XLAI,RGL,RSMIN,RSMAX . /1.,0.55,2.0,2.3,30.,100.,5000./ DATA TPREF,B1 /298.0,0.0016/ DATA C2 /60.0/ C ---------------------------------------------------------------------- C FOR LOW CROPS USE THESE VALUES FROM JACQUEMIN AND NOILHAN (90 BLM) C ---------------------------------------------------------------------- C DATA A1,A2,A3,A4,XLAI,RGL,RSMIN,RSMAX C . /1.,0.0055,0.55,2.0,2.3,100.,40.,5000./ tab = char(9) RCS = 0.0 RCT = 0.0 RCQ = 0.0 RCSOIL = 0.0 RC = 0.0 C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO INCOMING SOLAR RADIATION. C ---------------------------------------------------------------------- SOLAR = SOLARN/(1-ALBEDO) FF = A3*A4*SOLAR/(RGL*XLAI) RCS = (FF + RSMIN/RSMAX) / (A1 + FF) C rcs expression below from Pinty et al (89 JAM) C RCS = (A1 + A2*SOLAR) / ((RSMIN/RSMAX) + A2*SOLAR) C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO THE TEMPERATURE AT THE FIRST MODEL LEVEL. C ---------------------------------------------------------------------- RCT = 1.0 - B1*((TPREF-T2)**2.0) C rct expression below from Jarvis 76 and others C B4 = (THI-TOPT)/(TOPT-TLOW) C B3 = 1/((TOPT-TLOW)*(THI-TOPT)**B4) C RCT = B3*(T2-TLOW)*((THI-T2)**B4) C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT THE FIRST MODEL LEVEL. C ---------------------------------------------------------------------- P = SS(2)*PSFC CALL SVP(QS,ES,P,T2) RCQ = 1.0 - C2*(QS-W2) c WRITE (*,*) 'P,PSFC,SS(2)',p,psfc,ss(2) c WRITE (*,*) 'QS,ES,T2',QS,ES,T2 c WRITE (*,*) 'C2,W2,RCQ',C2,W2,RCQ RCQ = AMAX1(RCQ,0.25) c WRITE (*,*) 'MAX OF RCQ,0.25)',RCQ C rcq expression below from Noilhan and Planton (89 MWR) C P = SS(2)*PSFC C CALL SVP(QS,ES,P,T2) C EA = P*W2/(W2+0.622) C RCQ = 1-C1*(ES-EA) C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY. C DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP. C ---------------------------------------------------------------------- GX = (WSOIL(1) - TWILT) / (TSOREF - TWILT) IF (GX .GT. 1.) GX = 1. IF (GX .LT. 0.) GX = 0. PART(1) = (ZSOIL(1)/ZSOIL(NSOIL)) * GX DO 10 K = 2, NSOIL GX = (WSOIL(K) - TWILT) / (TSOREF - TWILT) IF (GX .GT. 1.) GX = 1. IF (GX .LT. 0.) GX = 0. PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NSOIL)) * GX 10 CONTINUE DO 20 K = 1, nsoil RCSOIL = RCSOIL+PART(K) 20 CONTINUE C ---------------------------------------------------------------------- C DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. C CONVERT CANOPY RESISTANCE (RC) TO PLANT COEFFICIENT (PC). C ---------------------------------------------------------------------- C RC = (RSMIN/XLAI)/(RCS*RCT*RCQ*RCSOIL) RC = RCMIN/(RCS*RCT*RCQ*RCSOIL) c PC = (RR + DELTA) / (RR * (1 + RC*CH) + DELTA) PC = (RR+SSK(1)-1.+DELTA) / . ((RR+SSK(1)-1.)*(1.+RC*CH)+DELTA) C ---------------------------------------------------------------------- C WRITE OUT EACH OF THE RC TERMS TO A FILE EVERY 10 TIME STEPS. C ---------------------------------------------------------------------- ircflg = ircflg+1 if (ircflg .eq. 10) then WRITE (12,1001) . TIMEIS,tab,rcs,tab,rct,tab,rcq,tab,rcsoil,tab,rc,tab,pc ircflg = 0 endif 1001 FORMAT (6(F10.4,A1),F10.4) RETURN END C ---------------------------------------------------------------------- C *** END SUBROUTINE CANRES *** C ----------------------------------------------------------------------