MODULE weather !- IPSL (2006) !- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC !--------------------------------------------------------------------- USE netcdf !- USE defprec USE ioipsl USE constantes USE parallel USE grid, ONLY : year,month,day,sec !- IMPLICIT NONE !- PRIVATE PUBLIC weathgen_main, weathgen_domain_size, weathgen_init, weathgen_read_file ! ! Only for root proc INTEGER, SAVE :: iim_file, jjm_file, llm_file, ttm_file INTEGER,DIMENSION(:,:),SAVE,ALLOCATABLE :: ncorr INTEGER,DIMENSION(:,:,:),SAVE,ALLOCATABLE :: icorr,jcorr INTEGER,SAVE :: i_cut, n_agg ! climatological wet days + anomaly (days/month) REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinwet ! climatological precipition + anomaly (mm/day) REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinprec ! climatological temp + anomaly (C) REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xint ! climatological relative humidity + anomaly (%) REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinq ! climatological wind speed + anomaly (m s-1) REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinwind ! climatological cloudiness + anomaly(%) REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xincld ! climatological temp range + anomaly(C) REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xintrng ! topography (m) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: xintopo ! latitudes of land points REAL,DIMENSION(:),SAVE,ALLOCATABLE :: lat_land !- ! daily values !- REAL,SAVE :: julian_last ! flag for wet day / dry day INTEGER,DIMENSION(:),SAVE,ALLOCATABLE :: iwet !- ! today's values (m0 means "minus 0") !- ! surface pressure (Pa) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: psurfm0 ! cloud fraction REAL,DIMENSION(:),SAVE,ALLOCATABLE :: cloudm0 ! maximum daily temperature (K) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tmaxm0 ! minimum daily temperature (K) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tminm0 ! daily average specific humidity (kg_h2o/kg_air) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: qdm0 ! daily average wind speed (m/sec) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: udm0 ! daily precitation (mm/day) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: precipm0 !- ! yesterday's values (m1 means "minus 1") !- ! surface pressure (Pa) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: psurfm1 ! cloud fraction REAL,DIMENSION(:),SAVE,ALLOCATABLE :: cloudm1 ! maximum daily temperature (K) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tmaxm1 ! minimum daily temperature (K) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tminm1 ! daily average specific humidity (kg_h2o/kg_air) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: qdm1 ! daily average wind speed (m/sec) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: udm1 ! daily precitation (mm/day) REAL,DIMENSION(:),SAVE,ALLOCATABLE :: precipm1 !- ! other !- ! statistical (0) or prescribed (1) daily values INTEGER,SAVE :: ipprec ! respect monthly precipitation LOGICAL,SAVE :: precip_exact INTEGER,DIMENSION(31,12),SAVE :: jour_precip ! max size of random seed INTEGER,PARAMETER :: seedsize_max = 300 LOGICAL,SAVE :: dump_weather CHARACTER(LEN=20),SAVE :: dump_weather_file LOGICAL,SAVE :: gathered INTEGER,SAVE :: dump_id ! ! Absolute zero REAL,PARAMETER :: zero_t=273.16 !- ! Parametres orbitaux: !- ! Eccentricity REAL,SAVE :: ecc ! Longitude of perihelie REAL,SAVE :: perh ! obliquity REAL,SAVE :: xob !- INTEGER,PARAMETER :: nmon = 12 ! CHARACTER(LEN=3),DIMENSION(12) :: cal = & & (/ 'JAN','FEB','MAR','APR','MAY','JUN', & & 'JUL','AUG','SEP','OCT','NOV','DEC' /) INTEGER,DIMENSION(12),SAVE :: ndaypm = & & (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) INTEGER,SAVE :: soldownid, rainfid, snowfid, lwradid, & & tairid, qairid, psolid, uid, vid, & & time_id, timestp_id !- ! Parameters for NETCDF : !- INTEGER,SAVE :: n_rtp = nf90_real4 !- ! Flag for dynamic allocations : INTEGER :: ALLOC_ERR !- ! Calendar type CHARACTER(LEN=20),SAVE :: calendar_str !- ! Land points index INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: kindex_w INTEGER, SAVE :: nbindex_w !- ! Plot of projection file => grid, number of column on terminal INTEGER, PARAMETER :: termcol = 100 !80 CONTAINS !- !=== !- SUBROUTINE daily & & (npoi, imonth, iday, cloud, tmax, tmin, precip, qd, ud, psurf) !--------------------------------------------------------------------- ! overview ! ! this routine generates daily weather conditions from monthly-mean ! climatic parameters ! ! specifically, this routine generates daily values of ! ! - daily total precipitation ! - daily maximum temperature ! - daily minimum temperature ! - daily average cloud cover ! - daily average relative humidity ! - daily average wind speed ! ! in order to generate daily weather conditions, the model uses ! a series of 'weather generator' approaches, ! which generate random combinations of weather conditions ! based upon the climatological conditions in general, ! this weather generator is based upon the so-called Richardson ! weather generator ! ! appropriate references include: ! ! Geng, S., F.W.T. Penning de Vries, and L. Supit, 1985: A simple ! method for generating rainfall data, Agricultural and Forest ! Meteorology, 36, 363-376. ! ! Richardson, C. W. and Wright, D. A., 1984: WGEN: A model for ! generating daily weather variables: U. S. Department of ! Agriculture, Agricultural Research Service. ! ! Richardson, C., 1981: Stochastic simulation of daily ! precipitation, temperature, and solar radiation. Water Resources ! Research 17, 182-190. !--------------------------------------------------------------------- !- ! in & out: global variables !- ! wet day / dry day flag ! INTEGER,INTENT(INOUT):: iwet(npoi) !- ! input !- ! total number of land points INTEGER,INTENT(IN) :: npoi INTEGER,INTENT(IN) :: imonth, iday !- ! output !- ! surface pressure (Pa) REAL,INTENT(OUT) :: psurf(npoi) ! cloud fraction REAL,INTENT(OUT) :: cloud(npoi) ! maximum daily temperature (K) REAL,INTENT(OUT) :: tmax(npoi) ! maximum daily temperature (K) REAL,INTENT(OUT) :: tmin(npoi) ! daily average specific humidity (kg_h2o/kg_air) REAL,INTENT(OUT) :: qd(npoi) ! daily average wind speed (m/sec) REAL,INTENT(OUT) :: ud(npoi) ! daily precitation (mm/day) REAL,INTENT(OUT) :: precip(npoi) !- ! local !- ! daily average temperature (K) REAL :: td(npoi) REAL,PARAMETER :: rair = 287. REAL,PARAMETER :: grav = 9.81 !!$ REAL,PARAMETER :: pi = 3.1415927 !- ! weather generator 'memory' matrix !- REAL,allocatable,save,dimension(:,:) :: xstore !- REAL :: ee(3), r(3), rr(npoi,3) REAL :: alpha(npoi), rndnum, pwd, pww REAL :: beta(npoi) REAL :: pwet(npoi) REAL :: rwork REAL :: omcloud, omqd, omtmax REAL :: cloudm, cloudw, cloudd REAL :: cloude(npoi), clouds(npoi) REAL :: tmaxd, tmaxw, tmaxm REAL :: tminm REAL :: tmins(npoi), tmaxs(npoi) REAL :: tmine(npoi), tmaxe(npoi) REAL :: qdm(npoi),qdd(npoi),qde(npoi),qdw(npoi),qdup(npoi),qdlow(npoi) INTEGER :: i,j,k REAL :: amn,b1,b2,b3,eud,rn,rn1,rn2,rn3,rn4,s1,s2,s12,x1,y, z(npoi) REAL :: aa(npoi),ab(npoi),tr1(npoi), tr2(npoi) REAL :: tdm,trngm,tdum REAL :: qsattd(npoi) INTEGER :: it1w, it2w REAL :: dt REAL :: rainpwd(npoi) INTEGER :: not_ok(npoi) INTEGER :: count_not_ok,count_not_ok_g LOGICAL,SAVE :: firstcall = .TRUE. INTEGER,save :: npoi0 REAL :: xx,e !- ! define autocorrelation matrices for Richardson generator ! ! note that this matrix should be based upon a statistical ! analysis of regional weather patterns ! ! for global simulations, we use 'nominal' values !- REAL, DIMENSION(3,3) :: a,b ! Warnings LOGICAL :: Warning_aa_ab(npoi), Warning_iwet(npoi) !--------------------------------------------------------------------- !- ! initial setup for daily climate calculations !- a(1,:) = (/ 0.600, 0.500, 0.005 /) a(2,:) = (/ 0.010, 0.250, 0.005 /) a(3,:) = (/ 0.020, 0.125, 0.250 /) !- b(1,:) = (/ 0.500, 0.250, -0.250 /) b(2,:) = (/ 0.000, 0.500, 0.250 /) b(3,:) = (/ 0.000, 0.000, 0.500 /) !- e = EXP(1.) ! GK240100 IF (firstcall) THEN firstcall = .FALSE. ALLOC_ERR=-1 ALLOCATE(xstore(npoi,3), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xstore : ",ALLOC_ERR STOP ENDIF xstore(:,:) = 0. npoi0 = npoi ELSE IF (npoi /= npoi0) THEN WRITE(numout,*) 'Domain size old, new: ',npoi0,npoi STOP 'WG Daily: Problem: Domain has changed since last call' ENDIF !- ! define working variables !- rwork = (grav/rair/0.0065) !- ! 'omega' parameters used to calculate differences in expected ! climatic parameters on wet and dry days ! ! following logic of weather generator used in the EPIC crop model ! ! omcloud -- cloud cover ! omqd -- humidity ! omtmax -- maximum temperature !- omcloud = 0.90 ! originally 0.90 omqd = 0.50 ! originally 0.50 omtmax = 0.75 ! originally 0.75 !- ! calculate weighting factors used in interpolating climatological ! monthly-mean input values to daily-mean values !- ! this is a simple linear interpolation technique that takes into ! account the length of each month !- IF (ipprec == 0) THEN IF (REAL(iday) < REAL(ndaypm(imonth)+1)/2.0) then it1w = imonth-1 it2w = imonth dt = (REAL(iday)-0.5)/ndaypm(imonth)+0.5 ELSE it1w = imonth it2w = imonth+1 dt = (REAL(iday)-0.5)/ndaypm(imonth)-0.5 ENDIF if (it1w < 1) it1w = 12 if (it2w > 12) it2w = 1 ELSE dt = -1. it1w = -1 it2w = -1 ENDIF !- IF (ipprec == 0) THEN !--- !-- use weather generator to create daily statistics !--- ! (1) determine if today will rain or not !--- !-- calculate monthly-average probability of rainy day !--- DO i=1,npoi pwet(i) = xinwet(i,imonth)/ndaypm(imonth) ENDDO !--- IF (.NOT.precip_exact) THEN !----- !---- (1.1) following Geng et al. !----- IF (is_root_prc) THEN CALL random_number (rndnum) ENDIF CALL bcast(rndnum) !----- DO i=1,npoi IF (xinprec(i,imonth) > 1.e-6) THEN !--------- !-------- implement simple first-order Markov-chain precipitation !-------- generator logic based on Geng et al. (1986), !-------- Richardson and Wright (1984), and Richardson (1981) !--------- !-------- basically, this allows for the probability of today being !-------- a wet day (a day with measureable precipitation) !-------- to be a function of what yesterday was (wet or dry) !--------- !-------- the logic here is that it is more likely that a wet day !-------- will follow another wet day -- allowing !-------- for 'storm events' to persist !--------- !-------- estimate the probability of a wet day after a dry day !--------- pwd = 0.75*pwet(i) !--------- !-------- estimate the probability of a wet day after a wet day !--------- pww = 0.25+pwd !--------- !-------- decide if today is a wet day or a dry day !-------- using a random number !--------- !-------- call random_number(rndnum) ! done before !--------- IF (iwet(i) == 0) then IF (rndnum <= pwd) iwet(i) = 1 ELSE IF (rndnum > pww) iwet(i) = 0 ENDIF ELSE iwet(i) = 0 ENDIF ENDDO ELSE !----- !---- (1.2) preserving the monthly number of precip days !---- and monthly precip !----- DO i=1,npoi IF (ABS(xinwet(i,imonth)) < 32.) THEN IF (xinprec(i,imonth) > 1.e-6) THEN IF ( jour_precip(iday,imonth) & & <= NINT(MAX(1.,xinwet(i,imonth))) ) THEN iwet(i) = 1 ELSE iwet(i) = 0 ENDIF ELSE iwet(i) = 0 ENDIF ENDIF ENDDO ENDIF !--- !-- (2) determine today's precipitation amount !--- IF (.not.precip_exact) THEN Warning_aa_ab(:)=.FALSE. Warning_iwet(:)=.FALSE. !----- !---- (2.1) following Geng et al. !----- aa(:) = zero ab(:) = zero tr2(:)= zero tr1(:)= zero beta(:) = un DO i=1,npoi !------- !------ initialize daily precipitation to zero !------- precip(i) = 0.0 !------- IF (xinprec(i,imonth) > 1.e-6) THEN !--------- !-------- if it is going to rain today !--------- IF (iwet(i) == 1) THEN !----------- !---------- calculate average rainfall per wet day !----------- rainpwd(i) = xinprec(i,imonth) & & *ndaypm(imonth)/MAX(0.1,xinwet(i,imonth)) !----------- !---------- randomly select a daily rainfall amount !---------- from a probability density function of rainfall !---------- !---------- method i -- !----------- !---------- use the following technique from Geng et al. and Richardson !---------- to distribute rainfall probabilities !----------- !---------- pick a random rainfall amount !---------- from a two-parameter gamma function distribution function !----------- !---------- estimate two parameters for gamma function !---------- (following Geng et al.) !----------- beta(i) = MAX(1.0,-2.16+1.83*rainpwd(i)) alpha(i) = rainpwd(i)/beta(i) !----------- !---------- determine daily precipitation amount !---------- from gamma distribution function !---------- (following WGEN code of Richardson and Wright (1984)) !----------- IF (ABS(1.-alpha(i)) < 1.e-6) THEN alpha(i) = 1.e-6*(alpha(i)/ABS(alpha(i))) ENDIF aa(i) = 1.0/alpha(i) ab(i) = 1.0/(1.0-alpha(i)) !----------- IF ( (ABS(aa(i)) < 1.e-6) .OR. (ABS(ab(i)) < 1.e-6) ) THEN Warning_aa_ab(:)=.TRUE. ENDIF tr1(i) = exp(-18.42/aa(i)) tr2(i) = exp(-18.42/ab(i)) ENDIF ELSE IF (iwet(i) == 1) THEN Warning_iwet(i)=.TRUE. ENDIF ENDIF ENDDO DO i=1,npoi IF ( Warning_aa_ab(i) ) THEN WRITE(numout,*) ' ATTENTION, aa ou ab:' WRITE(numout,*) ' aa, ab = ',aa(i),ab(i) WRITE(numout,*) ' alpha, rainpwd, beta =', & & alpha(i),rainpwd(i),beta(i) ENDIF IF ( Warning_iwet(i) ) THEN WRITE(numout,*) ' ATTENTION, iwet = 1 alors que xinprec = 0)' WRITE(numout,*) ' xinprec, iwet = ',xinprec(i,imonth),iwet(i) ENDIF ENDDO !----- where ( iwet(:) == 1 ) not_ok(:) = 1 elsewhere not_ok(:) = 0 endwhere !----- count_not_ok_g=0 count_not_ok=SUM(not_ok(:)) CALL reduce_sum(count_not_ok,count_not_ok_g) CALL bcast(count_not_ok_g) !- z(:)=0.0 DO WHILE (count_not_ok_g > 0) IF (is_root_prc) THEN CALL random_number (rn1) CALL random_number (rn2) ENDIF CALL bcast(rn1) CALL bcast(rn2) DO i=1,npoi IF ((iwet(i) == 1).AND.(not_ok(i) == 1)) then IF ( (rn1-tr1(i)) <= 0. ) THEN s1 = 0.0 ELSE s1 = rn1**aa(i) ENDIF !----------- IF ((rn2-tr2(i)) <= 0.) THEN s2 = 0.0 ELSE s2 = rn2**ab(i) ENDIF !----------- s12 = s1+s2 IF ((s12-1.0) <= 0.) THEN z(i) = s1/s12 not_ok(i) = 0 ENDIF ENDIF ENDDO count_not_ok_g=0 count_not_ok=SUM(not_ok(:)) CALL reduce_sum(count_not_ok,count_not_ok_g) CALL bcast(count_not_ok_g) ENDDO !----- IF (is_root_prc) THEN CALL random_number (rn3) ENDIF CALL bcast(rn3) ! WRITE(*,*) mpi_rank,"rn3=",rn3 !----- DO i=1,npoi IF (iwet(i) == 1) then precip(i) = -z(i)*log(rn3)*beta(i) ENDIF ENDDO !----- !---- method ii -- !----- !---- here we use a one-parameter Weibull distribution function !---- following the analysis of Selker and Haith (1990) !----- !---- Selker, J.S. and D.A. Haith, 1990: Development and testing !---- of single- parameter precipitation distributions, !---- Water Resources Research, 11, 2733-2740. !----- !---- this technique seems to have a significant advantage over other !---- means of generating rainfall distribution functions !----- !---- by calibrating the Weibull function to U.S. precipitation records, !---- Selker and Haith were able to establish the following relationship !----- !---- the cumulative probability of rainfall intensity x is given as: !----- !---- f(x) = 1.0-exp(-(1.191 x / rainpwd)**0.75) !----- !---- where x : rainfall intensity !---- rainpwd : rainfall per wet day !----- !---- using transformation method, take uniform deviate and convert !---- it to a random number weighted by the following Weibull function !----- !---- call random_number(rndnum) !----- !---- precip(i) = rainpwd / 1.191*(-log(1.0-rndnum))**1.333333 !----- !---- bound daily precipitation to "REAListic" range !----- DO i=1,npoi IF (iwet(i) == 1) THEN !--------- !-------- lower end is determined by definition of a 'wet day' !-------- (at least 0.25 mm of total precipitation) !--------- !-------- upper end is to prevent ibis from blowing up !--------- precip(i) = MAX(precip(i), 0.25) ! min = 0.25 mm/day precip(i) = MIN(precip(i),150.00) ! max = 150.00 mm/day ENDIF ENDDO ELSE !----- !---- (2.2) preserving the monthly number of precip days !---- and monthly precip !----- DO i=1,npoi !------- !------ Correction Nathalie. C'est abs(xinwet) < 32 qu'il faut tester !------ et non pas abs(xinprec(i,imonth)) < 32. !------- !------ IF ( (iwet(i) == 1).and.(abs(xinprec(i,imonth)) < 32.) ) THEN IF ( (iwet(i) == 1).and.(abs(xinwet(i,imonth)) < 32.) ) THEN precip(i) = xinprec(i,imonth)*REAL(ndaypm(imonth)) & & /REAL(NINT(MAX(1.,xinwet(i,imonth)))) ELSE precip(i) = 0. ENDIF ENDDO ENDIF !--- !-- (3) estimate expected minimum and maximum temperatures !--- DO i=1,npoi !----- !---- first determine the expected maximum and minimum temperatures !---- (from climatological means) for this day of the year !----- !---- mean daily mean temperature (K) tdm = xint(i,it1w)+dt*(xint(i,it2w)-xint(i,it1w))+zero_t !---- mean daily temperature range (K) trngm = xintrng(i,it1w)+dt*(xintrng(i,it2w)-xintrng(i,it1w)) !---- mean minimum and maximum temperatures tmaxm = tdm+0.56*trngm tminm = tdm-0.44*trngm !----- !---- modify maximum temperatures for wet and dry days !----- IF (pwet(i) /= 0.0) THEN tmaxd = tmaxm+pwet(i)*omtmax*trngm tmaxw = tmaxd- omtmax*trngm ELSE tmaxd = tmaxm tmaxw = tmaxm ENDIF !----- !---- set the 'expected' maximum and minimum temperatures for today !----- !---- note that the expected minimum temperatures are the same for !---- both wet and dry days !----- if (iwet(i) == 0) tmaxe(i) = tmaxd if (iwet(i) == 1) tmaxe(i) = tmaxw !----- tmine(i) = tminm !----- !---- estimate variability in minimum and maximum temperatures !----- !---- tmaxs : standard deviation in maximum temperature (K) !---- tmins : standard deviation in minimum temperature (K) !----- !---- Regression is based on analysis of 2-m air temperature data !---- from the NCEP/NCAR reanalysis (1958-1997) for 294 land points !---- over central North America !---- (24N-52N, 130W-60W, 0.5-degree resolution): !---- Daily max and min temperatures were calculated for each !---- land point from daily mean temperature and temperature range. !---- Anomalies were calculated !---- by subtracting similar max and min temperatures calculated from !---- monthly mean temperature and range (interpolated to daily). !---- The 40 years of anomalies were then binned by month !---- and the standard deviation calculated for each month. !---- The 294 x 12 standard deviations were then regressed !---- against the 3528 long-term monthly mean temperatures. !----- !---- note: the values are bound to be greater than 1.0 K !---- (at the very least they must be bound !---- so they don't go below zero) !----- tmaxs(i) = MAX(1.0,-0.0713*(tdm-zero_t)+4.89) tmins(i) = MAX(1.0,-0.1280*(tdm-zero_t)+5.73) ENDDO !--- !-- (4) estimate expected cloud cover !--- !--- !-- the formulation of dry and wet cloud cover has been !-- derived from the weather generator used in the epic crop model !--- DO i=1,npoi !----- !---- cloudm : mean cloud cover for today !---- cloudd : dry day cloud cover !---- cloudw : dry day cloud cover !---- cloude : expected cloud cover today !----- !---- mean cloud cover (%) !----- cloudm = xincld(i,it1w)+dt*(xincld(i,it2w)-xincld(i,it1w)) !----- !---- convert from percent to fraction !----- cloudm = cloudm/100.0 !----- !---- adjust cloud cover depending on dry day / rainy day !---- following logic of the EPIC weather generator code !----- IF (pwet(i) /= 0.0) THEN cloudd = (cloudm-pwet(i)*omcloud)/(1.0-pwet(i)*omcloud) cloudd = MIN(1.0,MAX(0.0,cloudd)) cloudw = (cloudm-(1.0-pwet(i))*cloudd)/pwet(i) ELSE cloudd = cloudm cloudw = cloudm ENDIF IF (iwet(i) == 0) cloude(i) = cloudd IF (iwet(i) == 1) cloude(i) = cloudw !----- !---- estimate variability in cloud cover for wet and dry days !---- following numbers proposed by Richardson !----- !---- clouds : standard deviation of cloud fraction !----- IF (iwet(i) == 0) clouds(i) = 0.24*cloude(i) IF (iwet(i) == 1) clouds(i) = 0.48*cloude(i) ENDDO !--- ! (5) determine today's temperatures and cloud cover using ! first-order serial autoregressive technique !--- !-- use the Richardson (1981) weather generator approach to simulate the !-- daily values of minimum / maximum temperature and cloud cover !--- !-- following the implementation of the Richardson WGEN weather !-- generator used in the EPIC crop model !--- !-- this approach uses a multivariate generator, which assumes that the !-- perturbation of minimum / maximum temperature and cloud cover are !-- normally distributed and that the serial correlation of each !-- variable may be described by a first-order autoregressive model !--- !-- generate standard deviates for weather generator !--- DO j=1,3 ee(j) = 9999. DO WHILE (ee(j) > 2.5) IF (is_root_prc) THEN CALL random_number (rn1) CALL random_number (rn2) ENDIF CALL bcast(rn1) CALL bcast(rn2) ee(j) = SQRT(-2.0*LOG(rn1))*COS(6.283185*rn2) ENDDO ENDDO !--- !-- zero out vectors !--- r(1:3) = 0.0 rr(1:npoi,1:3) = 0.0 !--- !-- update working vectors !--- do j=1,3 do k=1,3 r(j) = r(j)+b(j,k)*ee(j) enddo enddo !--- do j=1,3 do k=1,3 do i=1,npoi rr(i,j) = rr(i,j)+a(j,k)*xstore(i,k) enddo enddo enddo !--- !-- solve for x() perturbation vector and save current vector !-- into the xim1() storage vector (saved for each point) !--- do j=1,3 do i=1,npoi xstore(i,j) = r(j)+rr(i,j) enddo enddo !--- !-- determine today's minimum and maximum temperature !-- do i=1,npoi tmax(i) = tmaxe(i)+tmaxs(i)*xstore(i,1) tmin(i) = tmine(i)+tmins(i)*xstore(i,2) !----- !---- if tmin > tmax, then switch the two around !----- if (tmin(i) > tmax(i)) then tdum = tmax(i) tmax(i) = tmin(i) tmin(i) = tdum ENDIF !---- daily average temperature td(i) = 0.44*tmax(i)+0.56*tmin(i) !---- determine today's cloud cover cloud(i) = cloude(i)+clouds(i)*xstore(i,3) !---- constrain cloud cover to be between 0 and 100% cloud(i) = MAX(0.0,MIN(1.0,cloud(i))) enddo !--- !-- (6) estimate today's surface atmospheric pressure !--- do i=1,npoi !----- !---- simply a function of the daily average temperature !---- and topographic height -- nothing fancy here !----- psurf(i) = 101325.0*(td(i)/(td(i)+0.0065*xintopo(i)))**rwork enddo !--- !-- (7) estimate today's relative humidity !--- IF (is_root_prc) THEN CALL random_number (rn) ENDIF CALL bcast(rn) !--- CALL weathgen_qsat (npoi,td,psurf,qsattd) !--- do i=1,npoi !----- !---- the formulation of dry and wet relative humidities has been !---- derived from the weather generator used in the epic crop model !----- !---- qdm : mean relative humidity !---- qdd : dry day relative humidity !---- qdw : rainy day relative humidity !---- qde : expected relative humidity (based on wet/dry decision) !----- !---- mean relative humidity (%) qdm(i) = xinq(i,it1w)+dt*(xinq(i,it2w)-xinq(i,it1w)) !---- convert from percent to fraction qdm(i) = qdm(i)/100.0 !----- !---- adjust humidity depending on dry day / rainy day !---- following logic of the EPIC weather generator code !----- if (pwet(i) /= 0.0) then qdd(i) = (qdm(i)-pwet(i)*omqd)/(1.0-pwet(i)*omqd) if (qdd(i) < 0.2) then qdd(i) = 0.2 if (qdd(i) > qdm(i)) qdm(i) = qdd(i) ENDIF qdd(i) = MIN(1.0,qdd(i)) qdw(i) = (qdm(i)-(1.0-pwet(i))*qdd(i))/pwet(i) ELSE qdd(i) = qdm(i) qdw(i) = qdm(i) ENDIF !----- if (iwet(i) == 0) qde(i) = qdd(i) if (iwet(i) == 1) qde(i) = qdw(i) !----- !---- estimate lower and upper bounds of humidity distribution function !---- following logic of the EPIC weather generator code !----- xx = exp(qde(i)) qdup(i) = qde(i)+(1.0-qde(i))*xx/e qdlow(i) = qde(i)*(1.0-1./xx) !----- !---- randomlly select humidity from triangular distribution function !---- following logic of the EPIC weather generator code !----- !---- call random_number(rn) ! GK done before !----- y = 2.0/(qdup(i)-qdlow(i)) !----- b3 = qde(i)-qdlow(i) b2 = qdup(i)-qde(i) b1 = rn/y !----- x1 = y*b3/2.0 !----- if (rn > x1) then qd(i) = qdup(i)-sqrt (b2*b2-2.0*b2*(b1-0.5*b3)) ELSE qd(i) = qdlow(i)+sqrt (2.0*b1*b3) ENDIF !----- !---- adjust daily humidity to conserve monthly mean values !----- !---- note that this adjustment sometimes gives rise to humidity !---- values greater than 1.0 -- which is corrected below !----- amn = (qdup(i)+qde(i)+qdlow(i))/3.0 qd(i) = qd(i)*qde(i)/amn !----- !---- constrain daily average relative humidity !----- qd(i) = MAX(0.30,MIN(qd(i),0.99)) !----- !---- convert from relative humidity to specific humidity at !---- daily mean temperature !----- qd(i) = qd(i)*qsattd(i) enddo !--- !-- (8) estimate today's daily average wind speed !--- IF (is_root_prc) THEN CALL random_number (rn4) ENDIF CALL bcast(rn4) do i=1,npoi !----- !---- first estimate the expected daily average wind speed (from monthly !---- means) !----- eud = xinwind(i,it1w)+dt*(xinwind(i,it2w)-xinwind(i,it1w)) !----- !---- following logic of the EPIC weather generator !---- select random wind speed following this equation !----- !---- call random_number(rn4) !----- ud(i) = 1.13989*eud*(-log(rn4))**0.30 !---- constrain daily wind speeds to be between 2.5 and 10.0 m/sec ud(i) = MAX(2.5,MIN(ud(i),10.0)) ENDDO ELSE !--- !-- use REAL daily climate data !--- DO i=1,npoi !----- !---- use basic daily climate data, converting units !----- !---- daily total precipitation precip(i) = xinprec(i,imonth) !---- daily average temperatures td(i) = xint(i,imonth)+zero_t trngm = MIN(44.0,xintrng(i,imonth)) !----- tmax(i) = td(i)+0.56*trngm tmin(i) = td(i)-0.44*trngm !---- daily average cloud cover cloud(i) = xincld(i,imonth)*0.01 !---- daily average specific humidity qd(i) = xinq(i,imonth) !---- daily average wind speed ud(i) = xinwind(i,imonth) !----- !---- compute surface atmospheric pressure !----- psurf(i) = 101325.0*(td(i)/(td(i)+0.0065*xintopo(i)))**rwork ENDDO ENDIF !------------------- END SUBROUTINE daily !- !=== !- SUBROUTINE diurnal & & (npoi, nband, time, jday, plens, startp, endp, latitude, & & cloud, tmax, tmin, precip, qd, ud, psurf, & & fira, solad, solai, ua, ta, qa, raina, snowa, rh) !--------------------------------------------------------------------- IMPLICIT NONE !- ! input !- ! number of grid points INTEGER,INTENT(IN) :: npoi ! number of visible bands INTEGER,INTENT(IN) :: nband REAL,INTENT(IN) :: time INTEGER, INTENT(IN) :: jday REAL,INTENT(IN) :: plens,startp,endp ! latitude in degrees REAL,INTENT(IN) :: latitude(npoi) ! cloud fraction [0,1] REAL,INTENT(IN) :: cloud(npoi) ! maximum daily temperature (K) REAL,INTENT(IN) :: tmax(npoi) ! maximum daily temperature (K) REAL,INTENT(IN) :: tmin(npoi) ! daily precitation (mm/day) REAL,INTENT(IN) :: precip(npoi) ! daily average specific humidity (kg_h2o/kg_air) REAL,INTENT(IN) :: qd(npoi) ! daily average wind speed (m/sec) REAL,INTENT(IN) :: ud(npoi) ! surface pressure (Pa) REAL,INTENT(IN) :: psurf(npoi) !- ! output !- ! incoming ir flux (W m-2) REAL,INTENT(OUT) :: fira(npoi) ! direct downward solar flux (W m-2) REAL,INTENT(OUT) :: solad(npoi,nband) ! diffuse downward solar flux (W m-2) REAL,INTENT(OUT) :: solai(npoi,nband) ! wind speed (m s-1) REAL,INTENT(OUT) :: ua(npoi) ! air temperature (K) REAL,INTENT(OUT) :: ta(npoi) ! specific humidity (kg_h2o/kg_air) REAL,INTENT(OUT) :: qa(npoi) ! rainfall rate (mm/day) REAL,INTENT(OUT) :: raina(npoi) ! snowfall rate (mm/day) REAL,INTENT(OUT) :: snowa(npoi) ! relative humidity(%) REAL,INTENT(OUT) :: rh(npoi) !- ! local !- REAL,PARAMETER :: stef = 5.67051E-8 REAL,PARAMETER :: pi = 3.1415927 REAL,SAVE :: step REAL,PARAMETER :: pir = pi/180. REAL :: xl,so,xllp,xee,xse REAL :: xlam,dlamm,anm,ranm,ranv,anv,tls,rlam REAL :: sd,cd,deltar,delta,Dis_ST,ddt !- REAL :: coszen(npoi) ! cosine of solar zenith angle REAL :: rtime REAL :: orbit,angle,xdecl,xlat REAL :: sw,frac,gamma,qmin,qmax,qsa,emb,ea,ec,dtair,dtcloud,rn REAL :: trans(npoi), fdiffuse(npoi), qsatta(npoi), qsattmin(npoi) INTEGER :: i,ib INTEGER,SAVE :: npoi0 LOGICAL,SAVE :: firstcall = .TRUE. !--------------------------------------------------------------------- ! GK240100 IF (firstcall) THEN IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN step = 1.0 ELSE step = one_year/365.2425 ENDIF firstcall = .FALSE. npoi0 = npoi ELSE IF (npoi /= npoi0) THEN WRITE(numout,*) 'Domain size old, new: ',npoi0,npoi STOP 'WG Diurnal: Problem: Domain has changed since last call' ENDIF !- ! calendar and orbital calculations !- ! calculate time in hours rtime = time/3600.0 !- ! calculate the earth's orbital angle (around the sun) in radians orbit = 2.0*pi*REAL(jday)/365.2425 !- ! calculate the solar hour angle in radians angle = 2.0*pi*(rtime-12.0)/24.0 !- ! calculate the current solar declination angle ! ref: global physical climatology, hartmann, appendix a ! ! xdecl = 0.006918 & ! -0.399912*cos(orbit) & ! +0.070257*sin(orbit) & ! -0.006758*cos(2.0*orbit) & ! +0.000907*sin(2.0*orbit) & ! -0.002697*cos(3.0*orbit) & ! +0.001480*sin(3.0*orbit) ! ! calculate the effective solar constant, ! including effects of eccentricity ! ref: global physical climatology, hartmann, appendix a ! ! sw = 1370.*( 1.000110 & ! +0.034221*cos(orbit) & ! +0.001280*sin(orbit) & ! +0.000719*cos(2.0*orbit) & ! +0.000077*sin(2.0*orbit)) ! ! correction Marie-France Loutre ! ! orbital parameters ! ! ecc = 0.016724 ! perh = 102.04 ! xob = 23.446 !- xl = perh+180.0 ! so : sinus of obliquity so = sin(xob*pir) !- xllp = xl*pir xee = ecc*ecc xse = sqrt(1.0d0-xee) ! xlam : true long. sun for mean long. = 0 xlam = (ecc/2.0+ecc*xee/8.0d0)*(1.0+xse)*sin(xllp)-xee/4.0 & & *(0.5+xse)*sin(2.0*xllp)+ecc*xee/8.0*(1.0/3.0+xse) & & *sin(3.0*xllp) xlam =2.0d0*xlam/pir ! dlamm : mean long. sun for ma-ja dlamm =xlam+(jday-79)*step anm = dlamm-xl ranm = anm*pir xee = xee*ecc ! ranv : anomalie vraie (radian) ranv = ranm+(2.0*ecc-xee/4.0)*sin(ranm)+5.0/4.0*ecc*ecc & & *sin(2.0*ranm)+13.0/12.0*xee*sin(3.0*ranm) ! anv : anomalie vraie (degrees) anv = ranv/pir ! tls : longitude vraie (degrees) tls = anv+xl ! rlam : longitude vraie (radian) rlam = tls*pir ! sd and cd: cosinus and sinus of solar declination angle (delta) ! sinus delta = sin (obl)*sin(lambda) with lambda = real longitude ! (Phd. thesis of Marie-France Loutre, ASTR-UCL, Belgium, 1993) sd = so*sin(rlam) cd = sqrt(1.0d0-sd*sd) ! delta : Solar Declination (degrees, angle sun at equator) deltar = atan(sd/cd) delta = deltar/pir !- ! Eccentricity Effect !- Dis_ST =(1.0-ecc*ecc)/(1.0+ecc*cos(ranv)) ! ddt : 1 / normalized earth's sun distance ddt = 1.0/Dis_ST !- ! Insolation normal to the atmosphere (W/m2) !- sw = ddt *ddt *1370.d0 !- xdecl = deltar !- ! solar calculations !- do i=1,npoi !--- !-- calculate the latitude in radians !--- xlat = latitude(i)*pi/180.0 !--- !-- calculate the cosine of the solar zenith angle !--- coszen(i) = MAX(0.0, (sin(xlat)*sin(xdecl) & & + cos(xlat)*cos(xdecl)*cos(angle))) !--- !-- calculate the solar transmission through the atmosphere !-- using simple linear function of tranmission and cloud cover !--- !-- note that the 'cloud cover' data is typically obtained from !-- sunshine hours -- not direct cloud observations !--- !-- where, cloud cover = 1 - sunshine fraction !--- !-- different authors present different values for the slope and !-- intercept terms of this equation !--- !-- Friend, A: Parameterization of a global daily weather generator !-- for terrestrial ecosystem and biogeochemical modelling, !-- Ecological Modelling !--- !-- Spitters et al., 1986: Separating the diffuse and direct component !-- of global radiation and its implications for modeling canopy !-- photosynthesis, Part I: Components of incoming radiation, !-- Agricultural and Forest Meteorology, 38, 217-229. !--- !-- A. Friend : trans = 0.251+0.509*(1.0-cloud(i)) !-- Spitters et al. : trans = 0.200+0.560*(1.0-cloud(i)) !--- !-- we are using the values from A. Friend !--- trans(i) = 0.251+0.509*(1.0-cloud(i)) !--- !-- calculate the fraction of indirect (diffuse) solar radiation !-- based upon the cloud cover !--- !-- note that these relationships typically are measured for either !-- monthly or daily timescales, and may not be exactly appropriate !-- for hourly calculations -- however, in ibis, cloud cover is fixed !-- through the entire day so it may not make much difference !--- !-- method i -- !--- !-- we use a simple empirical relationships !-- from Nikolov and Zeller (1992) !--- !-- Nikolov, N. and K.F. Zeller, 1992: A solar radiation algorithm !-- for ecosystem dynamics models, Ecological Modelling, 61, 149-168. !--- fdiffuse(i) = 1.0045+trans(i)*( 0.0435+trans(i) & & *(-3.5227+trans(i)*2.6313)) !--- IF (trans(i) > 0.75) fdiffuse(i) = 0.166 !--- !-- method ii -- !--- !-- another method was suggested by Spitters et al. (1986) based on !-- long-term data from the Netherlands !-- !-- Spitters et al., 1986: Separating the diffuse and direct component !-- of global radiation and its implications for modeling canopy !-- photosynthesis, Part I: Components of incoming radiation, !-- Agricultural and Forest Meteorology, 38, 217-229. !--- !-- if ((trans == 0.00).and.(trans < 0.07)) then !-- fdiffuse = 1.0 !-- else if ((trans >= 0.07).and.(trans < 0.35)) then !-- fdiffuse = 1.0-2.3*(trans-0.07)**2 !-- else if ((trans >= 0.35).and.(trans < 0.75)) then !-- fdiffuse = 1.33-1.46*trans !-- ELSE !-- fdiffuse = 0.23 !-- endif !--- ENDDO !- ! do for each waveband !- DO ib=1,nband !--- !-- calculate the fraction in each waveband !--- !-- GK010200 IF (nband == 2) then frac = 0.46+0.08*REAL(ib-1) ELSE frac = 1./REAL(nband) ENDIF !--- DO i=1,npoi !----- !---- calculate the direct and indirect solar radiation !----- solad(i,ib) = sw*coszen(i)*frac*trans(i)*(1.-fdiffuse(i)) solai(i,ib) = sw*coszen(i)*frac*trans(i)*fdiffuse(i) ENDDO ENDDO !- ! temperature calculations !- !--- !-- assign hourly temperatures using tmax and tmin !-- following Environmental Biophysics, by Campbell and Norman, p.23 !--- !-- this function fits a fourier series to the diurnal temperature cycle !-- note that the maximum temperature occurs at 2:00 pm local solar time !--- !-- note that the daily mean value of gamma is 0.44, !-- so td = 0.44*tmax+0.56*tmin, instead of !-- td = 0.50*tmax+0.50*tmin !--- gamma = 0.44-0.46*SIN( pi/12.0*rtime+0.9) & +0.11*SIN(2.0*pi/12.0*rtime+0.9) DO i=1,npoi ta(i) = tmax(i)*gamma+tmin(i)*(1.0-gamma) ENDDO !- ! humidity calculations !- CALL weathgen_qsat (npoi,tmin,psurf,qsattmin) CALL weathgen_qsat (npoi,ta,psurf,qsatta) !- DO i=1,npoi !--- !-- adjust specific humidity against daily minimum temperatures !--- !-- To do this, qa is written as an approximate sine function !-- (same as ta) to preserve the daily mean specific humidity, !-- while also preventing rh from exceeding 99% at night !--- !-- Note that the daily mean RH is *not* preserved, and therefore the !-- output RH will be slightly different from what was read in. !--- !-- first adjust so that maximum RH cannot exceed 99% at night !--- qmin = MIN(qd(i),0.99*qsattmin(i)) qmax = (qd(i)-0.56*qmin)/0.44 !--- !-- if needed, adjust again to 99% at other times of the day (in which !-- case the daily mean *specific* humidity is also not preserved) !--- qsa = 0.99*qsatta(i) !--- !-- calculate the hourly specific humidity, using the above adjustments !--- qa(i) = MIN(qsa,qmax*gamma+qmin*(1.0-gamma)) !--- !-- calculate the hourly relative humidity !-- rh(i) = 100.*qa(i)/qsatta(i) ENDDO !- ! wind speed calculations !- IF (is_root_prc) THEN CALL random_number (rn) ENDIF CALL bcast(rn) !- DO i=1,npoi !--- !-- following logic of the EPIC weather generator !-- select random wind speed following this equation !--- !-- call random_number(rn) ! done before !--- ua(i) = 1.13989*ud(i)*(-log(rn))**0.30 !--- !-- fix wind speeds to always be above 2.5 m/sec and below 10.0 m/sec !--- ua(i) = MAX(2.5,MIN(10.0,ua(i))) ENDDO !- ! ir flux calculations !- DO i=1,npoi !--- !-- clear-sky emissivity as a function of water vapor pressure !-- and atmospheric temperature !--- !-- calculate the ir emissivity of the clear sky !-- using equation from idso (1981) water resources res., 17, 295-304 !--- emb = 0.01*(psurf(i)*qa(i)/(0.622+qa(i))) ea = 0.70+5.95e-5*emb*EXP(1500.0/ta(i)) !--- !-- assign the ir emissivity of clouds !-- (assume they are ~black in the ir) !--- ec = 0.950 !--- !-- assign the temperature difference of emissions (air+cloud) from !-- the surface air temperature !--- dtair = 0.0 dtcloud = 0.0 !--- !-- total downward ir is equal to the sum of: !--- !-- (1) clear sky contribution to downward ir radiation flux !-- (2) cloud contribution to downward ir radiation flux !--- fira(i) = (1.-cloud(i))*ea*stef*(ta(i)-dtair)**4 & & +cloud(i)*ec*stef*(ta(i)-dtcloud)**4 ENDDO !- ! snow and rain calculations !- DO i=1,npoi !--- !-- reset snow and rain to zero !--- snowa(i) = 0.0 raina(i) = 0.0 !--- !-- if precipitation event then calculate !--- IF (time >= startp .and. time < endp) then !----- !---- for rain / snow partitioning, make it all rain if !---- ta > 2.5 C and all snow if ta <= 2.5 C !----- !---- reference: !----- !---- Auer, A. H., 1974: The rain versus snow threshold temperatures, !---- Weatherwise, 27, 67. !----- IF (ta(i)-273.15 > 2.5) then raina(i) = precip(i)/plens ELSE snowa(i) = precip(i)/plens ENDIF ENDIF ENDDO !--------------------- END SUBROUTINE diurnal !- !=== !- SUBROUTINE weathgen_main & & (itau, istp, filename, force_id, iim, jjm, nband, & & rest_id, lrstread, lrstwrite, & & limit_west, limit_east, limit_north, limit_south, & & zonal_res, merid_res, lon, lat, date0, dt_force, & & kindex, nbindex, & & swdown, rainf, snowf, tair, u, v, qair, pb, lwdown) !--------------------------------------------------------------------- IMPLICIT NONE !- INTEGER,INTENT(IN) :: itau,istp CHARACTER(LEN=*),INTENT(IN) :: filename INTEGER,INTENT(IN) :: force_id INTEGER,INTENT(IN) :: iim,jjm INTEGER,INTENT(IN) :: nband INTEGER,INTENT(IN) :: rest_id LOGICAL,INTENT(IN) :: lrstread, lrstwrite REAL,INTENT(IN) :: limit_west,limit_east REAL,INTENT(IN) :: limit_north,limit_south REAL,INTENT(IN) :: zonal_res,merid_res REAL,DIMENSION(iim,jjm),INTENT(IN) :: lon,lat REAL,INTENT(IN) :: date0,dt_force !- INTEGER,DIMENSION(iim*jjm),INTENT(INOUT) :: kindex INTEGER,INTENT(INOUT) :: nbindex !- REAL,DIMENSION(iim,jjm),INTENT(OUT) :: & & tair,pb,qair,swdown,rainf,snowf, u,v,lwdown REAL, ALLOCATABLE, DIMENSION(:,:) :: & & tair_g,pb_g,qair_g,swdown_g,rainf_g,snowf_g, u_g,v_g,lwdown_g REAL,DIMENSION(nbindex) :: & & tairl,pbl,qairl,swdownl,rainfl,snowfl, ul,vl,lwdownl INTEGER :: i,j,ij !--------------------------------------------------------------------- IF (lrstread) THEN CALL weathgen_begin ( & & dt_force,itau, date0, & & rest_id,iim,jjm, & & lon,lat,tair,pb,qair,kindex,nbindex) RETURN ELSE CALL weathgen_get & & (itau, date0, dt_force, nbindex, nband, lat_land, & & swdownl, rainfl, snowfl, tairl, ul, vl, qairl, pbl, lwdownl) tair(:,:)=val_exp qair(:,:)=val_exp pb(:,:)=val_exp !!$ tair(:,:)=280. !!$ qair(:,:)=1.E-03 !!$ pb(:,:)=101325 DO ij=1,nbindex j = (((kindex(ij)-1)/iim) + 1) i = (kindex(ij) - (j-1)*iim) ! swdown(i,j)=swdownl(ij) rainf(i,j)=rainfl(ij) snowf(i,j)=snowfl(ij) tair(i,j)=tairl(ij) u(i,j)=ul(ij) v(i,j)=vl(ij) qair(i,j)=qairl(ij) pb(i,j)=pbl(ij) lwdown(i,j)=lwdownl(ij) ENDDO !--- IF (dump_weather) THEN ALLOCATE(tair_g(iim_g,jjm_g)) ALLOCATE(pb_g(iim_g,jjm_g)) ALLOCATE(qair_g(iim_g,jjm_g)) ALLOCATE(swdown_g(iim_g,jjm_g)) ALLOCATE(rainf_g(iim_g,jjm_g)) ALLOCATE(snowf_g(iim_g,jjm_g)) ALLOCATE(u_g(iim_g,jjm_g)) ALLOCATE(v_g(iim_g,jjm_g)) ALLOCATE(lwdown_g(iim_g,jjm_g)) CALL gather2D(tair, tair_g) CALL gather2D(pb, pb_g) CALL gather2D(qair, qair_g) CALL gather2D(swdown, swdown_g) CALL gather2D(rainf, rainf_g) CALL gather2D(snowf, snowf_g) CALL gather2D(u, u_g) CALL gather2D(v, v_g) CALL gather2D(lwdown, lwdown_g) IF (is_root_prc) THEN CALL weathgen_dump & & (itau, dt_force, iim_g, jjm_g, nbp_glo, index_g, lrstwrite, & & swdown_g, rainf_g, snowf_g, tair_g, u_g, v_g, qair_g, pb_g, lwdown_g) ENDIF ENDIF !--- IF (lrstwrite) THEN CALL weathgen_restwrite (rest_id,istp,iim,jjm,nbindex,kindex) ENDIF ENDIF !--------------------------- END SUBROUTINE weathgen_main !- !=== !- SUBROUTINE weathgen_init & & (filename,dt_force,force_id,iim,jjm, & & zonal_res,merid_res,lon,lat,kindex,nbindex,iind,jind) !--------------------------------------------------------------------- IMPLICIT NONE !- CHARACTER(LEN=*),INTENT(IN) :: filename REAL,INTENT(IN) :: dt_force INTEGER,INTENT(IN) :: force_id INTEGER,INTENT(IN) :: iim, jjm REAL,INTENT(IN) :: zonal_res,merid_res REAL,DIMENSION(iim,jjm),INTENT(IN) :: lon,lat !- INTEGER,DIMENSION(iim*jjm),INTENT(OUT) :: kindex INTEGER,INTENT(OUT) :: nbindex INTEGER,DIMENSION(iim),INTENT(OUT) :: iind INTEGER,DIMENSION(jjm),INTENT(OUT) :: jind !- REAL,PARAMETER :: fcrit = .5 REAL,DIMENSION(:),ALLOCATABLE :: lon_file, lon_temp REAL,DIMENSION(:,:),ALLOCATABLE :: nav_lon, nav_lat REAL,DIMENSION(:),ALLOCATABLE :: lat_file, lat_temp REAL,DIMENSION(:,:),ALLOCATABLE :: lsm_file REAL :: xextent_file, yextent_file, xres_file, yres_file INTEGER :: i, j, plotstep REAL,DIMENSION(iim,jjm) :: mask CHARACTER(LEN=1),DIMENSION(0:1) :: maskchar CHARACTER(LEN=30) :: var_name REAL :: x_cut REAL,DIMENSION(iim) :: tmp_lon REAL,DIMENSION(jjm) :: tmp_lat !--------------------------------------------------------------------- !- ! 0. messages, initialisations !- WRITE(numout,*) & & 'weathgen_init: Minimum land fraction on original grid =',fcrit CALL ioget_calendar(calendar_str) !- ! 1. on lit les longitudes et latitudes de la grille de depart ! ainsi que le masque terre-ocean !- CALL flinclo(force_id) CALL flininfo (filename,iim_file,jjm_file,llm_file,ttm_file,force_id) !- ALLOC_ERR=-1 ALLOCATE(nav_lon(iim_file,jjm_file), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of nav_lon : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(nav_lat(iim_file,jjm_file), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of nav_lat : ",ALLOC_ERR STOP ENDIF !- var_name='nav_lon' CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,nav_lon) var_name='nav_lat' CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,nav_lat) !- ALLOC_ERR=-1 ALLOCATE(lon_file(iim_file), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_file : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lat_file(jjm_file), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_file : ",ALLOC_ERR STOP ENDIF !- DO i=1,iim_file lon_file(i) = nav_lon(i,1) ENDDO DO j=1,jjm_file lat_file(j) = nav_lat(1,j) ENDDO !- ALLOC_ERR=-1 ALLOCATE(lsm_file(iim_file,jjm_file), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lsm_file : ",ALLOC_ERR STOP ENDIF !- var_name='lsmera' CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,lsm_file) !- ! 2. La resolution du modele ne doit pas etre superieure ! a celle de la grille de depart !- xextent_file = lon_file(iim_file)-lon_file(1) yextent_file = lat_file(1)-lat_file(jjm_file) xres_file = xextent_file/REAL(iim_file-1) yres_file = yextent_file/REAL(jjm_file-1) !- IF (xres_file > zonal_res) THEN WRITE(numout,*) 'Zonal resolution of model grid too fine.' WRITE(numout,*) 'Resolution of original data (deg): ', xres_file STOP ENDIF !- IF (yres_file > merid_res) THEN WRITE(numout,*) 'Meridional resolution of model grid too fine.' WRITE(numout,*) 'Resolution of original data (deg): ', yres_file STOP ENDIF !- ! 3. On verifie la coherence des coordonnees de depart et d'arrivee. ! Sinon, il faut deplacer une partie du monde (rien que ca). !- i_cut = 0 !- ALLOC_ERR=-1 ALLOCATE(lon_temp(iim_file), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_temp : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lat_temp(jjm_file), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_temp : ",ALLOC_ERR STOP ENDIF !- IF (lon(iim,1) > lon_file(iim_file)) THEN !-- A l'Est de la region d'interet, il n'y a pas de donnees !-- le bout a l'Ouest de la region d'interet est deplace de 360 deg !-- vers l'Est x_cut = lon(1,1) DO i=1,iim_file IF (lon_file(i) < x_cut) i_cut = i ENDDO IF ((i_cut < 1).OR.(i_cut >= iim_file)) THEN STOP 'Cannot find longitude for domain shift' ENDIF !--- lon_temp(1:iim_file-i_cut-1) = lon_file(i_cut:iim_file) lon_temp(iim_file-i_cut:iim_file) = lon_file(1:i_cut+1)+360. lon_file(:) = lon_temp(:) ELSEIF (lon(1,1) < lon_file(1)) THEN !-- A l'Ouest de la region d'interet, il n'y a pas de donnees !-- le bout a l'Est de la region d'interet est deplace de 360 deg !-- vers l'Ouest x_cut = lon(iim,1) DO i=1,iim_file IF (lon_file(i) < x_cut) i_cut = i ENDDO IF ( ( i_cut < 1 ) .OR. ( i_cut >= iim_file ) ) THEN STOP 'Cannot find longitude for domain shift' ENDIF !--- lon_temp(1:iim_file-i_cut-1) = lon_file(i_cut:iim_file) -360. lon_temp(iim_file-i_cut:iim_file) = lon_file(1:i_cut+1) lon_file(:) = lon_temp(:) ENDIF !- DEALLOCATE (lon_temp) DEALLOCATE (lat_temp) !- IF ( (lon(1,1) < lon_file(1)) & & .OR.( (lon(iim,1) > lon_file(iim_file)) & & .AND.(lon(iim,1) > lon_file(1)+360.) ) ) THEN WRITE(numout,*) lon(:,1) WRITE(numout,*) WRITE(numout,*) lon_file(:) STOP 'weathgen_init: cannot find coherent x-coordinates' ENDIF !- IF (i_cut /= 0) THEN CALL shift_field (iim_file,jjm_file,i_cut,lsm_file) ENDIF !- ! 4. Verification !- WRITE(numout,*) WRITE(numout,*) 'Input File: (Shifted) Global Land-Sea Mask' WRITE(numout,*) maskchar(0) = '-' maskchar(1) = 'X' plotstep = INT(REAL(iim_file-1)/termcol)+1 DO j=1,jjm_file,plotstep DO i=1,iim_file,plotstep WRITE(numout,'(a1,$)') maskchar(NINT(lsm_file(i,j))) ENDDO WRITE(numout,*) ENDDO WRITE(numout,*) !- ! 5. Prepare interpolation from fine grid land points to model grid !- ! 5.1 how many grid points of the original grid fall into one grid ! box of the model grid? !- n_agg = NINT((zonal_res/xres_file*merid_res/yres_file )+1.) !- ! au diable l'avarice ! !- n_agg = n_agg*2 !- ! 5.2 Allocate arrays where we store information about which ! (and how many) points on the original grid fall ! into which box on the model grid !- ALLOC_ERR=-1 ALLOCATE(ncorr(iim,jjm), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of ncorr : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(icorr(iim,jjm,n_agg), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of icorr : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(jcorr(iim,jjm,n_agg), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of jcorr : ",ALLOC_ERR STOP ENDIF !- ! 6. Land-Ocean Mask on model grid !- tmp_lon = lon(:,1) tmp_lat = lat(1,:) CALL mask_c_o & & (iim_file, jjm_file, lon_file, lat_file, lsm_file, fcrit, & & iim, jjm, zonal_res, merid_res, n_agg, tmp_lon, tmp_lat, & ! & iim, jjm, zonal_res, merid_res, n_agg, lon(:,1), lat(1,:), & & mask, ncorr, icorr, jcorr) !- WRITE(numout,*) 'Land-Sea Mask on Model Grid' WRITE(numout,*) plotstep = INT(REAL(iim-1)/termcol)+1 DO j=1,jjm,plotstep DO i=1,iim,plotstep WRITE(numout,'(a1,$)') maskchar(NINT(mask(i,j))) ENDDO WRITE(numout,*) ENDDO WRITE(numout,*) !- ! 7. kindex table. !- nbindex = 0 DO j=1,jjm DO i=1,iim IF (NINT(mask(i,j)) == 1) THEN nbindex = nbindex+1 kindex(nbindex) = (j-1)*iim+i ENDIF ENDDO ENDDO nbindex_w = nbindex ALLOC_ERR=-1 ALLOCATE(kindex_w(nbindex_w), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of kindex_w : ",ALLOC_ERR STOP ENDIF kindex_w(:)=kindex(1:nbindex) !- IF ( nbindex == 0 ) THEN WRITE(numout,*) 'Couillon! On est au plein milieu de l''ocean.' STOP 'Ou est-ce un bug?' ELSE WRITE(numout,*) 'Number of continental points: ',nbindex ENDIF !- ! 10. clean up !- DEALLOCATE (lon_file) DEALLOCATE (lat_file) DEALLOCATE (lsm_file) END SUBROUTINE weathgen_init SUBROUTINE weathgen_read_file & & (force_id,iim,jjm) !--------------------------------------------------------------------- IMPLICIT NONE !- INTEGER,INTENT(IN) :: force_id INTEGER,INTENT(IN) :: iim, jjm !- REAL,PARAMETER :: fcrit = .5 CHARACTER(LEN=30) :: var_name INTEGER,DIMENSION(iim*jjm) :: kindex INTEGER :: nbindex REAL,DIMENSION(iim*jjm) :: xchamp REAL,DIMENSION(iim*jjm,nmon) :: xchampm kindex(:)=0 #ifdef CPP_PARA nbindex=nbp_loc CALL scatter(index_g,kindex) kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g #else ! Copy saved land points index nbindex = nbindex_w kindex(1:nbindex_w) = kindex_w(:) #endif !- ALLOC_ERR=-1 ALLOCATE(lat_land(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_land : ",ALLOC_ERR STOP ENDIF !- ! 8 topography !- ALLOC_ERR=-1 ALLOCATE(xintopo(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xintopo : ",ALLOC_ERR STOP ENDIF !- var_name='altitude' CALL weather_read (force_id,var_name,iim_file,jjm_file,1,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchamp) xintopo(:)=xchamp(kindex(1:nbindex)) !- ! 9. Allocate and read the monthly fields !- ALLOC_ERR=-1 ALLOCATE(xinwet(nbindex,nmon), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xinwet : ",ALLOC_ERR STOP ENDIF var_name='prs' CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm) xinwet(:,:)=xchampm(kindex(1:nbindex),:) !- ALLOC_ERR=-1 ALLOCATE(xinprec(nbindex,nmon), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xinprec : ",ALLOC_ERR STOP ENDIF var_name='prm' CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm) xinprec(:,:)=xchampm(kindex(1:nbindex),:) !- ALLOC_ERR=-1 ALLOCATE(xint(nbindex,nmon), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xint : ",ALLOC_ERR STOP ENDIF var_name='t2m' CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm) xint(:,:)=xchampm(kindex(1:nbindex),:) !- ALLOC_ERR=-1 ALLOCATE(xinq(nbindex,nmon), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xinq : ",ALLOC_ERR STOP ENDIF var_name='r2m' CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm) xinq(:,:)=xchampm(kindex(1:nbindex),:) !- ALLOC_ERR=-1 ALLOCATE(xinwind(nbindex,nmon), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xinwind : ",ALLOC_ERR STOP ENDIF var_name='uv10m' CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm) xinwind(:,:)=xchampm(kindex(1:nbindex),:) !- ALLOC_ERR=-1 ALLOCATE(xincld(nbindex,nmon), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xincld : ",ALLOC_ERR STOP ENDIF var_name='tc' CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm) xincld(:,:)=xchampm(kindex(1:nbindex),:) !- ALLOC_ERR=-1 ALLOCATE(xintrng(nbindex,nmon), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xintrng : ",ALLOC_ERR STOP ENDIF var_name='t2a' CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, & iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm) xintrng(:,:)=xchampm(kindex(1:nbindex),:) !- ! 10. clean up !- IF (is_root_prc) THEN DEALLOCATE (ncorr) DEALLOCATE (icorr) DEALLOCATE (jcorr) ENDIF !- ! 12. Allocate space for daily mean values !- ALLOC_ERR=-1 ALLOCATE(iwet(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of iwet : ",ALLOC_ERR STOP ENDIF !- ALLOC_ERR=-1 ALLOCATE(psurfm0(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of psurfm0 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(cloudm0(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of cloudm0 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(tmaxm0(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of tmaxm0 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(tminm0(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of tminm0 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(qdm0(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of qdm0 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(udm0(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of udm0 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(precipm0(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of precipm0 : ",ALLOC_ERR STOP ENDIF !- ALLOC_ERR=-1 ALLOCATE(psurfm1(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of psurfm1 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(cloudm1(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of cloudm1 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(tmaxm1(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of tmaxm1 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(tminm1(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of tminm1 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(qdm1(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of qdm1 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(udm1(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of udm1 : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(precipm1(nbindex), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of precipm1 : ",ALLOC_ERR STOP ENDIF END SUBROUTINE weathgen_read_file SUBROUTINE weathgen_begin ( & & dt_force,itau, date0, & & rest_id,iim,jjm, & & lon,lat,tair,pb,qair,kindex,nbindex) !--------------------------------------------------------------------- IMPLICIT NONE !- REAL,INTENT(IN) :: dt_force, date0 INTEGER,INTENT(IN) :: itau INTEGER,INTENT(IN) :: rest_id INTEGER,INTENT(IN) :: iim, jjm REAL,DIMENSION(iim,jjm),INTENT(IN) :: lon,lat INTEGER,DIMENSION(iim*jjm),INTENT(OUT) :: kindex INTEGER,INTENT(OUT) :: nbindex !- REAL,DIMENSION(iim,jjm),INTENT(OUT) :: tair,pb,qair INTEGER :: i, j, ij REAL :: val_exp REAL,DIMENSION(iim*jjm) :: xchamp REAL,DIMENSION(iim_g*jjm_g) :: xchamp_g CHARACTER(LEN=30) :: var_name REAL,DIMENSION(1) :: jullasttab REAL,DIMENSION(seedsize_max) :: seed_in_file INTEGER,DIMENSION(:), ALLOCATABLE :: seed_in_proc INTEGER :: seedsize, iret INTEGER :: nlonid, nlatid, nlonid1, nlatid1, tdimid1, varid INTEGER :: ndim, dims(4) CHARACTER(LEN=30) :: assoc CHARACTER(LEN=70) :: str70 CHARACTER(LEN=80) :: stamp INTEGER :: yy_b, mm_b, dd_b, hh, mm REAL :: ss CHARACTER(LEN=10) :: today, att INTEGER :: nlandid1, nlandid, nlevid, nlevid1 REAL :: lev_max, lev_min REAL :: height_lev1 INTEGER :: imois REAL :: xx, td kindex(:)=0 #ifdef CPP_PARA nbindex=nbp_loc CALL scatter(index_g,kindex) kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g #else ! Copy saved land points index nbindex = nbindex_w kindex(1:nbindex_w) = kindex_w(:) #endif !- ! 13. Prescribed or statistical values? !- !Config Key = IPPREC !Config Desc = Use prescribed values !Config If = ALLOW_WEATHERGEN !Config Def = 0 !Config Help = If this is set to 1, the weather generator !Config uses the monthly mean values for daily means. !Config If it is set to 0, the weather generator !Config uses statistical relationships to derive daily !Config values from monthly means. ipprec = 0 CALL getin_p ('IPPREC',ipprec) WRITE(numout,*) 'IPPREC: ',ipprec !- ! 14. Do we want exact monthly precipitations even with ipprec=0 ? !- !Config Key = WEATHGEN_PRECIP_EXACT !Config Desc = Exact monthly precipitation !Config If = ALLOW_WEATHERGEN !Config Def = n !Config Help = If this is set to y, the weather generator !Config will generate pseudo-random precipitations !Config whose monthly mean is exactly the prescribed one. !Config In this case, the daily precipitation (for rainy !Config days) is constant (that is, some days have 0 precip, !Config the other days have precip=Precip_month/n_precip, !Config where n_precip is the prescribed number of rainy days !Config per month). precip_exact = .FALSE. CALL getin_p ('WEATHGEN_PRECIP_EXACT',precip_exact) WRITE(numout,*) 'PRECIP_EXACT: ',precip_exact !- IF (precip_exact) THEN !--- !-- preparer un tableau utilise pour determiner s'il pleut ou pas !-- (en fct. du nombre de jours de pluie par mois). !--- IF (is_root_prc) THEN DO imois=1,12 CALL permute (ndaypm(imois),jour_precip(:,imois)) ENDDO ENDIF CALL bcast(jour_precip) ENDIF !- ! Read Orbital Parameters !- !Config Key = ECCENTRICITY !Config Desc = Use prescribed values !Config If = ALLOW_WEATHERGEN !Config Def = 0.016724 ecc = 0.016724 CALL getin_p ('ECCENTRICITY',ecc) WRITE(numout,*) 'ECCENTRICITY: ',ecc ! !Config Key = PERIHELIE !Config Desc = Use prescribed values !Config If = ALLOW_WEATHERGEN !Config Def = 102.04 perh = 102.04 CALL getin_p ('PERIHELIE',perh) WRITE(numout,*) 'PERIHELIE: ',perh ! !Config Key = OBLIQUITY !Config Desc = Use prescribed values !Config If = ALLOW_WEATHERGEN !Config Def = 23.446 xob = 23.446 CALL getin_p ('OBLIQUITY',xob) WRITE(numout,*) 'OBLIQUITY: ',xob !- ! 15. Read restart file !- CALL ioget_expval (val_exp) !- var_name= 'julian' IF (is_root_prc) THEN CALL restget (rest_id,var_name,1,1,1,itau,.TRUE.,jullasttab) IF (jullasttab(1) == val_exp) THEN jullasttab(1) = itau2date(itau-1, date0, dt_force) ENDIF ENDIF CALL bcast(jullasttab) julian_last = jullasttab(1) !- var_name= 'seed' IF (is_root_prc) & CALL restget (rest_id,var_name,seedsize_max, & & 1,1,itau,.TRUE.,seed_in_file) CALL bcast(seed_in_file) IF (ALL(seed_in_file(:) == val_exp)) THEN !--- !-- there is no need to reinitialize the random number generator as !-- this does not seem to be a restart !--- CONTINUE ELSE !--- !-- reinitialize the random number generator !--- IF (is_root_prc) & CALL RANDOM_SEED( SIZE = seedsize ) CALL bcast(seedsize) IF (seedsize > seedsize_max) THEN STOP 'weathgen_begin: increase seedsize_max' ENDIF ALLOC_ERR=-1 ALLOCATE(seed_in_proc(seedsize), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of seed_in_proc : ",ALLOC_ERR STOP ENDIF seed_in_proc(1:seedsize) = NINT( seed_in_file(1:seedsize) ) CALL RANDOM_SEED (PUT = seed_in_proc) DEALLOCATE( seed_in_proc ) ENDIF !- var_name= 'iwet' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 0. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) iwet(:) = NINT(xchamp(kindex(1:nbindex))) !- var_name= 'psurfm0' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 100000. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) psurfm0(:) = xchamp(kindex(1:nbindex)) !- var_name= 'cloudm0' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = .5 ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) cloudm0(:) = xchamp(kindex(1:nbindex)) !- var_name= 'tmaxm0' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 285. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) tmaxm0(:) = xchamp(kindex(1:nbindex)) !- var_name= 'tminm0' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 275. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) tminm0(:) = xchamp(kindex(1:nbindex)) !- var_name= 'qdm0' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 1.E-03 ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) qdm0(:) = xchamp(kindex(1:nbindex)) !- var_name= 'udm0' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 2. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) udm0(:) = xchamp(kindex(1:nbindex)) !- var_name= 'precipm0' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 1. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) precipm0(:) = xchamp(kindex(1:nbindex)) !- var_name= 'psurfm1' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 100000. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) psurfm1(:) = xchamp(kindex(1:nbindex)) !- var_name= 'cloudm1' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = .5 ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) cloudm1(:) = xchamp(kindex(1:nbindex)) !- var_name= 'tmaxm1' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 285. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) tmaxm1(:) = xchamp(kindex(1:nbindex)) !- var_name= 'tminm1' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 275. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) tminm1(:) = xchamp(kindex(1:nbindex)) !- var_name= 'qdm1' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 1.E-03 ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) qdm1(:) = xchamp(kindex(1:nbindex)) !- var_name= 'udm1' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 2. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) udm1(:) = xchamp(kindex(1:nbindex)) !- var_name= 'precipm1' IF (is_root_prc) THEN CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g) IF (ALL(xchamp_g(:) == val_exp)) THEN xchamp_g(:) = 1. ENDIF ENDIF CALL scatter2D(xchamp_g,xchamp) precipm1(:) = xchamp(kindex(1:nbindex)) !- ! 16. We still need instantaneous tair, qair, and the surface pressure ! We take daily mean values read from the restart file !- !!$ tair(:,:)=280. !!$ qair(:,:)=1.E-03 !!$ pb(:,:)=101325 tair(:,:)=val_exp qair(:,:)=val_exp pb(:,:)=val_exp xx = 9.81/287./0.0065 DO ij=1,nbindex j = ((kindex(ij)-1)/iim) + 1 i = kindex(ij) - (j-1)*iim lat_land(ij) = lat(i,j) td = (tmaxm0(ij)+tminm0(ij))/2. tair(i,j) = td qair(i,j) = qdm1(ij) pb(i,j) = 101325.*(td/(td+0.0065*xintopo(ij)))**xx ENDDO !- ! 17. We can write a forcing file for Orchidee ! from this weather Generator run. !- !Config Key = DUMP_WEATHER !Config Desc = Write weather from generator into a forcing file !Config Def = n !Config Help = This flag makes the weather generator dump its ! generated weather into a forcing file which can ! then be used to get the same forcing on different ! machines. This only works correctly if there is ! a restart file (otherwise the forcing at the first ! time step is slightly wrong). dump_weather = .FALSE. CALL getin_p ('DUMP_WEATHER',dump_weather) !- IF (dump_weather .AND. is_root_prc) THEN !--- !-- Initialize the file !--- !Config Key = DUMP_WEATHER_FILE !Config Desc = Name of the file that contains ! the weather from generator !Config Def = 'weather_dump.nc' !Config If = DUMP_WEATHER !Config Help = dump_weather_file = 'weather_dump.nc' CALL getin ('DUMP_WEATHER_FILE',dump_weather_file) !--- !Config Key = DUMP_WEATHER_GATHERED !Config Desc = Dump weather data on gathered grid !Config Def = y !Config If = DUMP_WEATHER !Config Help = If 'y', the weather data are gathered ! for all land points. gathered = .TRUE. CALL getin ('DUMP_WEATHER_GATHERED',gathered) !--- iret = NF90_CREATE (TRIM(dump_weather_file),NF90_CLOBBER,dump_id) !--- !-- Dimensions !--- iret = NF90_DEF_DIM (dump_id,'x',iim_g,nlonid1) iret = NF90_DEF_DIM (dump_id,'y',jjm_g,nlatid1) iret = NF90_DEF_DIM (dump_id,'z', 1,nlevid1) !--- IF (gathered) THEN iret = NF90_DEF_DIM (dump_id,'land',nbp_glo,nlandid1) ENDIF iret = NF90_DEF_DIM (dump_id,'tstep',NF90_UNLIMITED,tdimid1) !--- !-- Coordinate variables !--- dims(1:2) = (/ nlonid1, nlatid1 /) !--- iret = NF90_DEF_VAR (dump_id,'nav_lon',n_rtp,dims(1:2),nlonid) iret = NF90_PUT_ATT (dump_id,nlonid,'units',"degrees_east") iret = NF90_PUT_ATT (dump_id,nlonid,'valid_min',MINVAL(lon_g)) iret = NF90_PUT_ATT (dump_id,nlonid,'valid_max',MAXVAL(lon_g)) iret = NF90_PUT_ATT (dump_id,nlonid,'long_name',"Longitude") !--- iret = NF90_DEF_VAR (dump_id,'nav_lat',n_rtp,dims(1:2),nlatid) iret = NF90_PUT_ATT (dump_id,nlatid,'units',"degrees_north") iret = NF90_PUT_ATT (dump_id,nlatid,'valid_min',MINVAL(lat_g)) iret = NF90_PUT_ATT (dump_id,nlatid,'valid_max',MAXVAL(lat_g)) iret = NF90_PUT_ATT (dump_id,nlatid,'long_name',"Latitude") !--- height_lev1 = 10. CALL getin ('HEIGHT_LEV1',height_lev1) lev_min = height_lev1 lev_max = height_lev1 !--- iret = NF90_DEF_VAR (dump_id,'level',n_rtp,(/ nlevid1 /),nlevid) iret = NF90_PUT_ATT (dump_id,nlevid,'units',"m") iret = NF90_PUT_ATT (dump_id,nlevid,'valid_min',lev_min) iret = NF90_PUT_ATT (dump_id,nlevid,'valid_max',lev_max) iret = NF90_PUT_ATT (dump_id,nlevid,'long_name',"Vertical levels") !--- IF (gathered) THEN iret = NF90_DEF_VAR (dump_id,'land',NF90_INT,(/ nlandid1 /),nlandid) iret = NF90_PUT_ATT (dump_id,nlandid,'compress',"y x") ENDIF !--- !-- Store the time axes. !--- iret = NF90_DEF_VAR (dump_id,'time',n_rtp,tdimid1,time_id) yy_b=0 mm_b=1 dd_b=1 hh=00 mm=00 ss=0. WRITE (str70,7000) yy_b, mm_b, dd_b, hh, mm, INT(ss) iret = NF90_PUT_ATT (dump_id,time_id,'units',TRIM(str70)) iret = NF90_PUT_ATT (dump_id,time_id,'calendar',TRIM(calendar_str)) iret = NF90_PUT_ATT (dump_id,time_id,'title','Time') iret = NF90_PUT_ATT (dump_id,time_id,'long_name','Time axis') WRITE(str70,7001) yy_b, cal(mm_b), dd_b, hh, mm, INT(ss) iret = NF90_PUT_ATT (dump_id,time_id,'time_origin',TRIM(str70)) !--- !-- Time steps !--- iret = NF90_DEF_VAR (dump_id,'timestp',NF90_INT,tdimid1,timestp_id) WRITE(str70,7002) yy_b, mm_b, dd_b, hh, mm, INT(ss) iret = NF90_PUT_ATT (dump_id,timestp_id,'units',TRIM(str70)) iret = NF90_PUT_ATT (dump_id,timestp_id,'title','Time steps') iret = NF90_PUT_ATT (dump_id,timestp_id,'tstep_sec',dt_force) iret = NF90_PUT_ATT & & (dump_id,timestp_id,'long_name','Time step axis') WRITE(str70,7001) yy_b, cal(mm_b), dd_b, hh, mm, INT(ss) iret = NF90_PUT_ATT (dump_id,timestp_id,'time_origin',TRIM(str70)) !--- 7000 FORMAT('seconds since ',I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 7001 FORMAT(' ',I4.4,'-',A3,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 7002 FORMAT('timesteps since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) !--- !-- The variables in the file are declared and their attributes !-- written. !--- IF (gathered) THEN ndim = 2 dims(1:2) = (/ nlandid1, tdimid1 /) assoc = 'time (nav_lat nav_lon)' ELSE ndim = 3 dims(1:3) = (/ nlonid1, nlatid1, tdimid1 /) assoc = 'time nav_lat nav_lon' ENDIF !--- iret = NF90_DEF_VAR (dump_id,'SWdown',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX') iret = NF90_PUT_ATT (dump_id,varid,'units','W/m^2') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Surface incident shortwave radiation') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) soldownid = varid !--- iret = NF90_DEF_VAR (dump_id,'Rainf',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX') iret = NF90_PUT_ATT (dump_id,varid,'units','kg/m^2s') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Rainfall rate') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) rainfid = varid !--- iret = NF90_DEF_VAR (dump_id,'Snowf',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX') iret = NF90_PUT_ATT (dump_id,varid,'units','kg/m^2s') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Snowfall rate') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) snowfid = varid !--- iret = NF90_DEF_VAR (dump_id,'LWdown',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX') iret = NF90_PUT_ATT (dump_id,varid,'units','W/m^2') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Surface incident longwave radiation') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) lwradid = varid !--- iret = NF90_DEF_VAR (dump_id,'PSurf',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX') iret = NF90_PUT_ATT (dump_id,varid,'units','Pa') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Surface pressure') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) psolid = varid !--- !-- 3D Variables to be written !--- IF (gathered) THEN ndim = 3 dims(1:3) = (/ nlandid1, nlevid1, tdimid1 /) assoc = 'time level (nav_lat nav_lon)' ELSE ndim = 4 dims(1:4) = (/ nlonid1, nlatid1, nlevid1, tdimid1 /) assoc = 'time level nav_lat nav_lon' ENDIF !--- iret = NF90_DEF_VAR (dump_id,'Tair',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX') iret = NF90_PUT_ATT (dump_id,varid,'units','K') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Near surface air temperature') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) tairid = varid !--- iret = NF90_DEF_VAR (dump_id,'Qair',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX') iret = NF90_PUT_ATT (dump_id,varid,'units','kg/kg') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Near surface specific humidity') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) qairid = varid !--- iret = NF90_DEF_VAR (dump_id,'Wind_N',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX') iret = NF90_PUT_ATT (dump_id,varid,'units','m/s') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Near surface northward wind component') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) uid = varid !--- iret = NF90_DEF_VAR (dump_id,'Wind_E',n_rtp,dims(1:ndim),varid) iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX') iret = NF90_PUT_ATT (dump_id,varid,'units','m/s') iret = NF90_PUT_ATT (dump_id,varid,'long_name', & & 'Near surface eastward wind component') iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc)) iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba) vid = varid !--- !-- Global attributes !--- CALL DATE_AND_TIME (today, att) stamp = "WG, date: "//TRIM(today)//" at "//TRIM(att) iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'Conventions',"GDT 1.2") iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'file_name', & & TRIM(dump_weather_file)) iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'production',TRIM(stamp)) !--- !-- Finish the definition phase !--- iret = NF90_ENDDEF (dump_id) !--- !-- Write coordinates !--- iret = NF90_PUT_VAR (dump_id,nlonid,lon_g) IF (iret /= NF90_NOERR) THEN WRITE(numout,*) iret CALL ipslerr (3,'weathgen_begin', & & 'Could not put variable nav_lon in the file : ', & & TRIM(dump_weather_file),'(Solution ?)') ENDIF iret = NF90_PUT_VAR (dump_id,nlatid,lat_g) IF (iret /= NF90_NOERR) THEN WRITE(numout,*) iret CALL ipslerr (3,'weathgen_begin', & & 'Could not put variable nav_lat in the file : ', & & TRIM(dump_weather_file),'(Solution ?)') ENDIF iret = NF90_PUT_VAR (dump_id,nlevid,height_lev1) IF (iret /= NF90_NOERR) THEN WRITE(numout,*) iret CALL ipslerr (3,'weathgen_begin', & & 'Could not put variable level in the file : ', & & TRIM(dump_weather_file),'(Solution ?)') ENDIF !--- IF (gathered) THEN iret = NF90_PUT_VAR (dump_id,nlandid,index_g(1:nbp_glo)) IF (iret /= NF90_NOERR) THEN WRITE(numout,*) iret CALL ipslerr (3,'weathgen_begin', & & 'Could not put variable land in the file : ', & & TRIM(dump_weather_file),'(Solution ?)') ENDIF ENDIF !--- ENDIF ! dump_weather !----------------------------- END SUBROUTINE weathgen_begin !- !=== !- SUBROUTINE weathgen_get & & (itau, date0, dt_force, nbindex, nband, lat, & & swdown, raina, snowa, tair, u, v, qair, psurf, lwdown) !--------------------------------------------------------------------- IMPLICIT NONE ! number of time step INTEGER,INTENT(IN) :: itau ! date when itau was 0 REAL,INTENT(IN) :: date0 ! time step (s) REAL,INTENT(IN) :: dt_force ! number of land points INTEGER,INTENT(IN) :: nbindex ! number of visible bands INTEGER,INTENT(IN) :: nband ! latitude (deg) REAL,DIMENSION(nbindex),INTENT(IN) :: lat !- REAL,DIMENSION(nbindex,nband),INTENT(OUT) :: swdown REAL,DIMENSION(nbindex),INTENT(OUT) :: raina, snowa REAL,DIMENSION(nbindex),INTENT(OUT) :: tair REAL,DIMENSION(nbindex),INTENT(OUT) :: u,v REAL,DIMENSION(nbindex),INTENT(OUT) :: qair REAL,DIMENSION(nbindex),INTENT(OUT) :: psurf REAL,DIMENSION(nbindex),INTENT(OUT) :: lwdown !- REAL,DIMENSION(nbindex) :: cloud, tmax, tmin, precipd, qd, ud REAL,DIMENSION(nbindex) :: rh REAL,DIMENSION(nbindex,nband) :: solai, solad REAL :: julian, jur REAL :: x INTEGER :: yy, mm, dd REAL :: ss, plens, time !--------------------------------------------------------------------- !- ! 1. get a reduced julian day !- julian = itau2date(itau-1, date0, dt_force) !SZ, test: solar noon at 12 o'clock! ! julian = itau2date(itau, date0, dt_force) CALL ju2ymds (julian, yy, mm, dd, ss) CALL ymds2ju (yy,1,1,0.0, jur) julian = julian-jur CALL ju2ymds (julian, yy, mm, dd, ss) !- ! 2. daily values !- IF (INT(julian_last) /= INT(julian)) THEN !-- today's values become yesterday's values cloudm1(:) = cloudm0(:) tmaxm1(:) = tmaxm0(:) tminm1(:) = tminm0(:) precipm1(:) = precipm0(:) qdm1(:) = qdm0(:) udm1(:) = udm0(:) psurfm1(:) = psurfm0(:) !-- we have to get new daily values !!$ WRITE(*,*) mpi_rank, "weathgen_get : date ",yy, mm, dd, ss !!$ WRITE(*,*) mpi_rank, "weathgen_get : grid date ",year, month, day, sec CALL daily (nbindex, mm, dd, cloudm0, tmaxm0, tminm0, & & precipm0, qdm0, udm0, psurfm0) ENDIF !- ! 3. interpolate daily values ! (otherwise we get ugly temperature jumps at midnight) !- x = (julian-INT(julian)) !- cloud(:) = (1.-x)*cloudm1(:)+x*cloudm0(:) tmax(:) = (1.-x)*tmaxm1(:)+x*tmaxm0(:) tmin(:) = (1.-x)*tminm1(:)+x*tminm0(:) precipd(:) = (1.-x)*precipm1(:)+x*precipm0(:) qd(:) = (1.-x)*qdm1(:)+x*qdm0(:) ud(:) = (1.-x)*udm1(:)+x*udm0(:) psurf(:) = (1.-x)*psurfm1(:)+x*psurfm0(:) !- ! 4. read instantaneous values !- plens = one_day/dt_force time = (julian-REAL(INT(julian)))*one_day !- CALL diurnal (nbindex, nband, time, NINT(julian), plens, 0., one_day, & & lat, cloud, tmax, tmin, precipd, qd, ud, psurf, & & lwdown, solad, solai, u, tair, qair, raina, snowa, rh) !- raina(:) = raina(:)/dt_force snowa(:) = snowa(:)/dt_force !- swdown(:,:) = solad(:,:)+solai(:,:) !- v(:) = 0. !- ! 5. Store date !- julian_last = julian !-------------------------- END SUBROUTINE weathgen_get !- !=== !- SUBROUTINE weathgen_restwrite (rest_id,itau,iim,jjm,nbindex,kindex) !--------------------------------------------------------------------- IMPLICIT NONE !- INTEGER,INTENT(IN) :: rest_id,itau,iim,jjm,nbindex INTEGER,DIMENSION(iim*jjm),INTENT(IN) :: kindex !- CHARACTER(LEN=30) :: var_name INTEGER :: i,j,ij REAL,DIMENSION(1) :: jullasttab REAL,DIMENSION(seedsize_max) :: seed_in_file INTEGER,DIMENSION(:),ALLOCATABLE :: seed_in_proc INTEGER :: seedsize REAL :: rndnum REAL,DIMENSION(iim*jjm) :: xchamp REAL,DIMENSION(iim_g*jjm_g) :: xchamp_g !--------------------------------------------------------------------- var_name= 'julian' jullasttab(:) = julian_last IF (is_root_prc) CALL restput (rest_id, var_name, 1, 1, 1, itau, jullasttab) !- IF (is_root_prc) THEN CALL RANDOM_SEED( SIZE = seedsize ) IF (seedsize > seedsize_max) THEN STOP 'weathgen_restwrite: increase seedsize_max' ENDIF ENDIF CALL bcast(seedsize) IF (is_root_prc) THEN ALLOC_ERR=-1 ALLOCATE(seed_in_proc(seedsize), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of seed_in_proc : ",ALLOC_ERR STOP ENDIF !- CALL RANDOM_SEED (GET = seed_in_proc) !- seed_in_file(1:seedsize) = REAL(seed_in_proc(1:seedsize)) !- ! fill in the seed up to seedsize_max ! (useful in the case we restart on ! a machine with a longer seed vector: ! we do not want a degenerated seed) !- DO i=seedsize+1,seedsize_max CALL RANDOM_NUMBER( rndnum ) seed_in_file(i) = 100000.*rndnum ENDDO ENDIF CALL bcast (seed_in_file) !- IF (is_root_prc) THEN DEALLOCATE( seed_in_proc ) !- var_name= 'seed' CALL restput (rest_id,var_name,seedsize_max,1,1,itau,seed_in_file) ENDIF !- xchamp(:) = val_exp !!$ DO j=1,jjm !!$ DO i=1,iim !!$ ij = (j-1)*iim+i !!$ xchamp(i,j) = REAL(iwet(ij)) !!$ ENDDO !!$ ENDDO DO ij=1,nbindex xchamp(kindex(ij)) = REAL(iwet(ij)) ENDDO var_name= 'iwet' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = psurfm0(ij) ENDDO var_name= 'psurfm0' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = cloudm0(ij) ENDDO var_name= 'cloudm0' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = tmaxm0(ij) ENDDO var_name= 'tmaxm0' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = tminm0(ij) ENDDO var_name= 'tminm0' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = qdm0(ij) ENDDO var_name= 'qdm0' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = udm0(ij) ENDDO var_name= 'udm0' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = precipm0(ij) ENDDO var_name= 'precipm0' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = psurfm1(ij) ENDDO var_name= 'psurfm1' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = cloudm1(ij) ENDDO var_name= 'cloudm1' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = tmaxm1(ij) ENDDO var_name= 'tmaxm1' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = tminm1(ij) ENDDO var_name= 'tminm1' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = qdm1(ij) ENDDO var_name= 'qdm1' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = udm1(ij) ENDDO var_name= 'udm1' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !- DO ij=1,nbindex xchamp(kindex(ij)) = precipm1(ij) ENDDO var_name= 'precipm1' CALL gather2D(xchamp,xchamp_g) IF (is_root_prc) THEN CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g) ENDIF !-------------------------------- END SUBROUTINE weathgen_restwrite !- !=== !- SUBROUTINE weather_read & & (force_id,nomvar,iim_file,jjm_file,n3,i_cut, & & iim,jjm,n_agg,ncorr,icorr,jcorr,champout) !--------------------------------------------------------------------- IMPLICIT NONE !- INTEGER,INTENT(IN) :: force_id CHARACTER(LEN=*),INTENT(IN) :: nomvar INTEGER,INTENT(IN) :: iim_file,jjm_file INTEGER,INTENT(IN) :: n3 INTEGER,INTENT(IN) :: i_cut INTEGER,INTENT(IN) :: iim,jjm INTEGER,INTENT(IN) :: n_agg INTEGER,DIMENSION(:,:),INTENT(IN) :: ncorr INTEGER,DIMENSION(:,:,:),INTENT(IN) :: icorr,jcorr !- REAL,DIMENSION(iim*jjm,n3),INTENT(OUT) :: champout !- REAL,DIMENSION(iim_file,jjm_file,n3) :: champ_file REAL,ALLOCATABLE,DIMENSION(:,:) :: champout_g INTEGER :: i,j,ij,l,m !--------------------------------------------------------------------- WRITE(numout,*) 'Lecture ',TRIM(nomvar) !- IF (is_root_prc) THEN ALLOCATE(champout_g(iim_g*jjm_g,n3)) IF ( n3 == 1 ) THEN CALL flinget (force_id,nomvar(1:LEN_TRIM(nomvar)), & & iim_file, jjm_file, 0, 0, 1, 1, champ_file) ELSE DO l=1,n3 CALL flinget & & (force_id,nomvar(1:LEN_TRIM(nomvar)), & & iim_file, jjm_file, 0, n3, l, l, champ_file(:,:,l)) ENDDO ENDIF ! shift if necessary IF (i_cut /= 0) THEN DO l=1,n3 CALL shift_field (iim_file,jjm_file,i_cut,champ_file(:,:,l)) ENDDO ENDIF ! interpolate onto the model grid DO l=1,n3 DO j=1,jjm_g DO i=1,iim_g ij = i+(j-1)*iim_g champout_g(ij,l) = 0. DO m=1,ncorr(i,j) champout_g(ij,l) = champout_g(ij,l) & & +champ_file(icorr(i,j,m),jcorr(i,j,m),l) ENDDO champout_g(ij,l) = champout_g(ij,l)/REAL(ncorr(i,j)) ENDDO ENDDO ENDDO !!$ DO l=1,n3 !!$ DO j=1,jjm_g !!$ WRITE(numout,*) j,(/ ( champout_g((j-1)*iim_g+i,l), i=1,iim_g ) /) !!$ ENDDO !!$ ENDDO ELSE ALLOCATE(champout_g(1,1)) ENDIF !!$ CALL scatter2D(champout_g,champout) #ifndef CPP_PARA champout(:,:)=champout_g(:,:) #else CALL scatter2D_rgen(champout_g,champout,n3) #endif !!$ DO l=1,n3 !!$ DO j=1,jjm !!$ WRITE(numout,*) j,(/ ( champout((j-1)*iim_g+i,l), i=1,iim_g ) /) !!$ ENDDO !!$ ENDDO !---------------------------- END SUBROUTINE weather_read !- !=== !- SUBROUTINE weathgen_dump & & (itau, dt_force, iim, jjm, nbindex, kindex, lrstwrite, & & swdown, rainf, snowf, tair, u, v, qair, pb, lwdown ) !--------------------------------------------------------------------- IMPLICIT NONE !- INTEGER,INTENT(IN) :: itau REAL,INTENT(IN) :: dt_force INTEGER,INTENT(IN) :: iim,jjm INTEGER,INTENT(IN) :: nbindex INTEGER,DIMENSION(iim*jjm),INTENT(IN) :: kindex LOGICAL,INTENT(IN) :: lrstwrite REAL,DIMENSION(iim*jjm),INTENT(IN) :: & & swdown, rainf, snowf, tair, u, v, qair, pb, lwdown !- INTEGER :: iret,ndim INTEGER,DIMENSION(4) :: corner,edges REAL,DIMENSION(iim*jjm) :: var_gather !--------------------------------------------------------------------- !- ! time dimension !- iret = NF90_PUT_VAR (dump_id,timestp_id,(/ REAL(itau) /), & & start=(/ itau /),count=(/ 1 /)) iret = NF90_PUT_VAR (dump_id,time_id,(/ REAL(itau)*dt_force /), & & start=(/ itau /),count=(/ 1 /)) !- ! 2D variables: pas de dimension verticale !- IF (gathered) THEN ndim = 2 corner(1:2) = (/ 1, itau /) edges(1:2) = (/ nbindex, 1 /) ELSE ndim = 3 corner(1:3) = (/ 1, 1, itau /) edges(1:3) = (/ iim, jjm, 1 /) ENDIF !- CALL gather_weather (iim*jjm,nbindex,kindex,swdown, var_gather) iret = NF90_PUT_VAR (dump_id,soldownid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) CALL gather_weather (iim*jjm,nbindex,kindex,rainf, var_gather) iret = NF90_PUT_VAR (dump_id,rainfid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) CALL gather_weather (iim*jjm,nbindex,kindex,snowf, var_gather) iret = NF90_PUT_VAR (dump_id,snowfid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) CALL gather_weather (iim*jjm,nbindex,kindex,pb, var_gather) iret = NF90_PUT_VAR (dump_id,psolid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) CALL gather_weather (iim*jjm,nbindex,kindex,lwdown, var_gather) iret = NF90_PUT_VAR (dump_id,lwradid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) !- ! 3D variables !- IF (gathered) THEN ndim = 3 corner(1:3) = (/ 1, 1, itau /) edges(1:3) = (/ nbindex, 1, 1 /) ELSE ndim = 4 corner(1:4) = (/ 1, 1, 1, itau /) edges(1:4) = (/ iim, jjm, 1, 1 /) ENDIF !- CALL gather_weather (iim*jjm,nbindex,kindex,u, var_gather) iret = NF90_PUT_VAR (dump_id,uid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) CALL gather_weather (iim*jjm,nbindex,kindex,v, var_gather) iret = NF90_PUT_VAR (dump_id,vid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) CALL gather_weather (iim*jjm,nbindex,kindex,tair, var_gather) iret = NF90_PUT_VAR (dump_id,tairid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) CALL gather_weather (iim*jjm,nbindex,kindex,qair, var_gather) iret = NF90_PUT_VAR (dump_id,qairid, var_gather, & & start=corner(1:ndim), count=edges(1:ndim)) !- IF (lrstwrite) THEN iret = NF90_CLOSE (dump_id) ENDIF !--------------------------- END SUBROUTINE weathgen_dump !- !=== !- SUBROUTINE gather_weather (iimjjm, nbindex, kindex, var_in, var_out) !--------------------------------------------------------------------- IMPLICIT NONE !- INTEGER,INTENT(IN) :: iimjjm,nbindex INTEGER,DIMENSION(iimjjm),INTENT(IN) :: kindex REAL,DIMENSION(iimjjm),INTENT(IN) :: var_in !- REAL,DIMENSION(iimjjm),INTENT(OUT) :: var_out !- INTEGER :: i LOGICAL,SAVE :: firstcall = .TRUE. INTEGER,SAVE :: nb_outside INTEGER,ALLOCATABLE,SAVE,DIMENSION(:) :: outside !--------------------------------------------------------------------- IF (firstcall) THEN !--- !-- determine which points are not in the computational domain and !-- create a mask for these points !--- firstcall = .FALSE. ALLOC_ERR=-1 ALLOCATE(outside(iimjjm), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of outside : ",ALLOC_ERR STOP ENDIF outside(:) = 0. nb_outside = 0 DO i=1,iimjjm IF ( ALL( kindex(:) /= i ) ) THEN nb_outside = nb_outside+1 outside(nb_outside) = i ENDIF ENDDO ENDIF !- IF ( gathered ) THEN DO i=1,nbindex var_out(i) = var_in(kindex(i)) ENDDO ELSE var_out(:) = var_in(:) DO i=1,nb_outside var_out(outside(i)) = undef_sechiba ENDDO ENDIF !-------------------- END SUBROUTINE gather_weather !- !=== !- SUBROUTINE shift_field (im,jm,i_cut,champ) !--------------------------------------------------------------------- INTEGER,INTENT(IN) :: im,jm,i_cut REAL,DIMENSION(im,jm),INTENT(INOUT) :: champ !- REAL,DIMENSION(im,jm) :: champ_temp !--------------------------------------------------------------------- IF ( (i_cut >= 1).AND.(i_cut <= im) ) THEN champ_temp(1:im-i_cut-1,:) = champ(i_cut:im,:) champ_temp(im-i_cut:im,:) = champ(1:i_cut+1,:) champ(:,:) = champ_temp(:,:) ENDIF !------------------------- END SUBROUTINE shift_field !- !=== !- SUBROUTINE weathgen_domain_size & & (limit_west,limit_east,limit_north,limit_south, & & zonal_res,merid_res,iim,jjm) !--------------------------------------------------------------------- IMPLICIT NONE !- REAL,INTENT(INOUT) :: limit_west,limit_east,limit_north,limit_south REAL,INTENT(IN) :: zonal_res,merid_res INTEGER,INTENT(OUT) :: iim,jjm !--------------------------------------------------------------------- IF (limit_west > limit_east) limit_east = limit_east+360. !- IF ( (limit_west >= limit_east) & & .OR.(limit_east > 360.) & & .OR.(limit_west < -180.) & & .OR.(limit_east-limit_west > 360.) ) THEN WRITE(numout,*) 'PROBLEME longitudes.' WRITE(numout,*) 'Limites Ouest, Est: ',limit_west,limit_east STOP ENDIF !- IF ( (limit_south < -90.) & & .OR.(limit_north > 90.) & & .OR.(limit_south >= limit_north ) ) THEN WRITE(numout,*) 'PROBLEME latitudes.' WRITE(numout,*) 'Limites Nord, Sud: ',limit_north,limit_south STOP ENDIF !- IF ( (zonal_res <= 0. ) & & .OR.(zonal_res > limit_east-limit_west) ) THEN WRITE(numout,*) 'PROBLEME resolution zonale.' WRITE(numout,*) 'Limites Ouest, Est, Resolution: ', & & limit_west,limit_east,zonal_res STOP ENDIF !- IF ( (merid_res <= 0.) & & .OR.(merid_res > limit_north-limit_south) ) THEN WRITE(numout,*) 'PROBLEME resolution meridionale.' WRITE(numout,*) 'Limites Nord, Sud, Resolution: ', & & limit_north,limit_south,merid_res STOP ENDIF !- iim = NINT(MAX((limit_east-limit_west)/zonal_res,1.)) jjm = NINT(MAX((limit_north-limit_south)/merid_res,1.)) !- WRITE(numout,*) 'Domain size: iim, jjm = ', iim, jjm !---------------------------------- END SUBROUTINE weathgen_domain_size !- !=== !- FUNCTION tsatl (t) RESULT (tsat) !--------------------------------------------------------------------- ! statement functions tsatl,tsati are used below so that lowe's ! polyomial for liquid is used if t gt 273.16, or for ice if ! t lt 273.16. also impose range of validity for lowe's polys. !--------------------------------------------------------------------- REAL,INTENT(IN) :: t REAL :: tsat !--------------------------------------------------------------------- tsat = MIN(100.,MAX(t-zero_t,0.)) !----------------- END FUNCTION tsatl !- !=== !- FUNCTION tsati (t) RESULT (tsat) !--------------------------------------------------------------------- ! statement functions tsatl,tsati are used below so that lowe's ! polyomial for liquid is used if t gt 273.16, or for ice if ! t lt 273.16. also impose range of validity for lowe's polys. !--------------------------------------------------------------------- REAL,INTENT(IN) :: t REAL :: tsat !--------------------------------------------------------------------- tsat = MAX(-60.,MIN(t-zero_t,0.)) !----------------- END FUNCTION tsati !- !=== !- FUNCTION esat (t) RESULT (esatout) !--------------------------------------------------------------------- ! statement function esat is svp in n/m**2, with t in deg k. ! (100 * lowe's poly since 1 mb = 100 n/m**2.) !--------------------------------------------------------------------- REAL,INTENT(IN) :: t REAL :: esatout REAL :: x !- ! polynomials for svp(t), d(svp)/dt over water and ice are from ! lowe(1977),jam,16,101-103. !- REAL,PARAMETER :: & & asat0 = 6.1078000, asat1 = 4.4365185e-1, asat2 = 1.4289458e-2, & & asat3 = 2.6506485e-4, asat4 = 3.0312404e-6, asat5 = 2.0340809e-8, & & asat6 = 6.1368209e-11, & & bsat0 = 6.1091780, bsat1 = 5.0346990e-1, bsat2 = 1.8860134e-2, & & bsat3 = 4.1762237e-4, bsat4 = 5.8247203e-6, bsat5 = 4.8388032e-8, & & bsat6 = 1.8388269e-10 !--------------------------------------------------------------------- IF (t >= zero_t) THEN x = asat0 ELSE x = bsat0 ENDIF !- esatout = 100.* & ( x & +tsatl(t)*(asat1+tsatl(t)*(asat2+tsatl(t)*(asat3 & +tsatl(t)*(asat4+tsatl(t)*(asat5+tsatl(t)* asat6))))) & +tsati(t)*(bsat1+tsati(t)*(bsat2+tsati(t)*(bsat3 & +tsati(t)*(bsat4+tsati(t)*(bsat5+tsati(t)* bsat6))))) ) !---------------- END FUNCTION esat !- !=== !- FUNCTION qsat (e,p) RESULT (qsatout) !--------------------------------------------------------------------- ! statement function qsat is saturation specific humidity, ! with svp e and ambient pressure p in n/m**2. impose an upper ! limit of 1 to avoid spurious values for very high svp ! and/or small p !--------------------------------------------------------------------- REAL, INTENT(IN) :: e,p REAL :: qsatout !--------------------------------------------------------------------- qsatout = 0.622*e/MAX(p-(1.0-0.622)*e,0.622*e) !---------------- END FUNCTION qsat !- !=== !- SUBROUTINE weathgen_qsat (npoi,t,p,qsat) !--------------------------------------------------------------------- ! vectorized version of functions esat and qsat. ! statement function esat is svp in n/m**2, with t in deg k. ! (100 * lowe's poly since 1 mb = 100 n/m**2.) !--------------------------------------------------------------------- INTEGER,INTENT(IN) :: npoi REAL,DIMENSION(npoi),INTENT(IN) :: t,p REAL,DIMENSION(npoi),INTENT(OUT):: qsat !- REAL,DIMENSION(npoi) :: x, tl, ti, e !- ! polynomials for svp(t), d(svp)/dt over water and ice ! are from lowe(1977),jam,16,101-103. !- REAL,PARAMETER :: & & asat0 = 6.1078000, asat1 = 4.4365185e-1, asat2 = 1.4289458e-2, & & asat3 = 2.6506485e-4, asat4 = 3.0312404e-6, asat5 = 2.0340809e-8, & & asat6 = 6.1368209e-11, & & bsat0 = 6.1091780, bsat1 = 5.0346990e-1, bsat2 = 1.8860134e-2, & & bsat3 = 4.1762237e-4, bsat4 = 5.8247203e-6, bsat5 = 4.8388032e-8, & & bsat6 = 1.8388269e-10 !--------------------------------------------------------------------- WHERE (t(:) > zero_t) x(:) = asat0 ELSEWHERE x(:) = bsat0 ENDWHERE !- tl(:) = MIN(100.,MAX(t(:)-zero_t,0.)) ti(:) = MAX(-60.,MIN(t(:)-zero_t,0.)) !- e(:) = 100.* & ( x(:) & +tl(:)*(asat1+tl(:)*(asat2+tl(:)*(asat3 & +tl(:)*(asat4+tl(:)*(asat5+tl(:)* asat6))))) & +ti(:)*(bsat1+ti(:)*(bsat2+ti(:)*(bsat3 & +ti(:)*(bsat4+ti(:)*(bsat5+ti(:)* bsat6))))) ) !- qsat(:) = 0.622*e(:)/MAX(p(:)-(1.0-0.622)*e(:),0.622*e(:)) !--------------------------- END SUBROUTINE weathgen_qsat !- !=== !- SUBROUTINE mask_c_o & & (imdep, jmdep, xdata, ydata, mask_in, fcrit, & & imar, jmar, zonal_res, merid_res, n_agg, x, y, mask, & & ncorr, icorr, jcorr) !--------------------------------------------------------------------- ! z.x.li (le 01/04/1994) : ! A partir du champ de masque, on fabrique ! un champ indicateur (masque) terre/ocean ! terre:1; ocean:0 !--------------------------------------------------------------------- INTEGER :: imdep,jmdep REAL :: xdata(imdep),ydata(jmdep) REAL :: mask_in(imdep,jmdep) REAL :: fcrit INTEGER :: imar,jmar REAL :: zonal_res,merid_res INTEGER :: n_agg REAL :: x(imar),y(jmar) REAL, INTENT(OUT) :: mask(imar,jmar) INTEGER :: ncorr(imar,jmar) INTEGER,DIMENSION(imar,jmar,n_agg) :: icorr,jcorr !- INTEGER i, j, ii, jj REAL a(imar),b(imar),c(jmar),d(jmar) INTEGER num_tot(imar,jmar), num_oce(imar,jmar) REAL,ALLOCATABLE :: distans(:) INTEGER ij_proche(1),i_proche,j_proche !- INTEGER,DIMENSION(imar,jmar) :: ncorr_oce , ncorr_land INTEGER,DIMENSION(imar,jmar,n_agg) :: & & icorr_oce, jcorr_oce , icorr_land, jcorr_land !- INTEGER imdepp REAL,ALLOCATABLE :: xdatap(:) REAL,ALLOCATABLE :: mask_inp(:,:) LOGICAL :: extend !--------------------------------------------------------------------- ncorr(:,:) = 0 icorr(:,:,:) = -1; jcorr(:,:,:) = -1 ncorr_land(:,:) = 0 icorr_land(:,:,:) = -1; jcorr_land(:,:,:) = -1 ncorr_oce(:,:) = 0 icorr_oce(:,:,:) = -1; jcorr_oce(:,:,:) = -1 ! do we have to extend the domain (-x...-x+360)? IF ( xdata(1)+360.-xdata(imdep) > 0.01 ) THEN extend = .TRUE. imdepp = imdep+1 ELSE extend = .FALSE. imdepp = imdep ENDIF !- ALLOC_ERR=-1 ALLOCATE(xdatap(imdepp), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of xdatap : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(mask_inp(imdepp,jmdep), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of mask_inp : ",ALLOC_ERR STOP ENDIF !- xdatap(1:imdep) = xdata(1:imdep) mask_inp(1:imdep,:) = mask_in(1:imdep,:) !- IF (extend) THEN xdatap(imdepp) = xdatap(1)+360. mask_inp(imdepp,:) = mask_inp(1,:) ENDIF !- ALLOC_ERR=-1 ALLOCATE(distans(imdepp*jmdep), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of distans : ",ALLOC_ERR STOP ENDIF ! Definition des limites des boites de la grille d'arrivee. IF (imar > 1) THEN a(1) = x(1)-(x(2)-x(1))/2.0 b(1) = (x(1)+x(2))/2.0 DO i=2,imar-1 a(i) = b(i-1) b(i) = (x(i)+x(i+1))/2.0 ENDDO a(imar) = b(imar-1) b(imar) = x(imar)+(x(imar)-x(imar-1))/2.0 ELSE a(1) = x(1)-zonal_res/2. b(1) = x(1)+zonal_res/2. ENDIF !- IF (jmar > 1) THEN c(1) = y(1)-(y(2)-y(1))/2.0 d(1) = (y(1)+y(2))/2.0 DO j=2,jmar-1 c(j) = d(j-1) d(j) = (y(j)+y(j+1))/2.0 ENDDO c(jmar) = d(jmar-1) d(jmar) = y(jmar)+(y(jmar)-y(jmar-1))/2.0 ELSE c(1) = y(1)-merid_res/2. d(1) = y(1)+merid_res/2. ENDIF !- num_oce(1:imar,1:jmar) = 0 num_tot(1:imar,1:jmar) = 0 !- ! ..... Modif P. Le Van ( 23/08/95 ) .... !- DO ii=1,imar DO jj=1,jmar DO i=1,imdepp IF ( ( (xdatap(i)-a(ii) >= 1.e-5) & & .AND.(xdatap(i)-b(ii) <= 1.e-5) ) & & .OR.( (xdatap(i)-a(ii) <= 1.e-5) & & .AND.(xdatap(i)-b(ii) >= 1.e-5) ) ) THEN DO j=1,jmdep IF ( ( (ydata(j)-c(jj) >= 1.e-5) & & .AND.(ydata(j)-d(jj) <= 1.e-5) ) & & .OR.( (ydata(j)-c(jj) <= 1.e-5) & & .AND.(ydata(j)-d(jj) >= 1.e-5) ) ) THEN num_tot(ii,jj) = num_tot(ii,jj)+1 IF (mask_inp(i,j) < 0.5) THEN num_oce(ii,jj) = num_oce(ii,jj)+1 !-------------- on a trouve un point oceanique. On le memorise. ncorr_oce(ii,jj) = ncorr_oce(ii,jj)+1 IF ((i == imdepp).AND.extend) THEN icorr_oce(ii,jj,ncorr_oce(ii,jj)) = 1 ELSE icorr_oce(ii,jj,ncorr_oce(ii,jj)) = i ENDIF jcorr_oce(ii,jj,ncorr_oce(ii,jj)) = j ELSE !-------------- on a trouve un point continental. On le memorise. ncorr_land(ii,jj) = ncorr_land(ii,jj)+1 IF ((i == imdepp).AND.extend) THEN icorr_land(ii,jj,ncorr_land(ii,jj)) = 1 ELSE icorr_land(ii,jj,ncorr_land(ii,jj)) = i ENDIF jcorr_land(ii,jj,ncorr_land(ii,jj)) = j ENDIF ENDIF ENDDO ENDIF ENDDO ENDDO ENDDO !- DO i=1,imar DO j=1,jmar IF (num_tot(i,j) > 0) THEN IF ( ( (num_oce(i,j) == 0) & & .AND.(num_tot(i,j) > 0) ) & & .OR.( (num_oce(i,j) > 0) & & .AND.( REAL(num_oce(i,j)) & & <= REAL(num_tot(i,j))*(1.-fcrit) ) ) ) THEN mask(i,j) = 1. ncorr(i,j) = ncorr_land(i,j) icorr(i,j,:) = icorr_land(i,j,:) jcorr(i,j,:) = jcorr_land(i,j,:) ELSE mask(i,j) = 0. ncorr(i,j) = ncorr_oce(i,j) icorr(i,j,:) = icorr_oce(i,j,:) jcorr(i,j,:) = jcorr_oce(i,j,:) ENDIF ELSE CALL dist_sphe(x(i),y(j),xdatap,ydata,imdepp,jmdep,distans) ij_proche(:) = MINLOC(distans) j_proche = (ij_proche(1)-1)/imdepp+1 i_proche = ij_proche(1)-(j_proche-1)*imdepp mask(i,j) = mask_inp(i_proche,j_proche) IF ( (i_proche == imdepp).AND.extend) i_proche=1 ncorr(i,j) = 1 icorr(i,j,1) = i_proche jcorr(i,j,1) = j_proche ENDIF ENDDO ENDDO !---------------------- END SUBROUTINE mask_c_o !- !=== !- SUBROUTINE dist_sphe (rf_lon,rf_lat,rlon,rlat,im,jm,distance) !--------------------------------------------------------------------- ! Auteur: Laurent Li (le 30/12/1996) ! ! Ce programme calcule la distance minimale (selon le grand cercle) ! entre deux points sur la terre !--------------------------------------------------------------------- !- ! Input: !- INTEGER im, jm ! dimensions REAL rf_lon ! longitude du point de reference (degres) REAL rf_lat ! latitude du point de reference (degres) REAL rlon(im), rlat(jm) ! longitude et latitude des points !- ! Output: !- REAL distance(im,jm) ! distances en metre !- REAL rlon1, rlat1 REAL rlon2, rlat2 REAL dist REAL pa, pb, p, pi INTEGER i,j !- REAL radius PARAMETER (radius=6371229.) !--------------------------------------------------------------------- pi = 4.0*ATAN(1.0) !- DO j=1,jm DO i=1,im rlon1=rf_lon rlat1=rf_lat rlon2=rlon(i) rlat2=rlat(j) pa = pi/2.0-rlat1*pi/180.0 ! dist. entre pole n et point a pb = pi/2.0-rlat2*pi/180.0 ! dist. entre pole n et point b !----- p = (rlon1-rlon2)*pi/180.0 ! angle entre a et b (leurs meridiens) !----- dist = ACOS( COS(pa)*COS(pb)+SIN(pa)*SIN(pb)*COS(p)) dist = radius*dist distance(i,j) = dist ENDDO ENDDO !----------------------- END SUBROUTINE dist_sphe !- !=== !- SUBROUTINE permute (n,ordre) !--------------------------------------------------------------------- INTEGER,INTENT(IN) :: n INTEGER,DIMENSION(n),INTENT(OUT) :: ordre !- INTEGER,DIMENSION(n) :: restant INTEGER :: ipique, i, n_rest REAL :: rndnum !--------------------------------------------------------------------- n_rest = n restant(:) = (/ (i, i=1,n) /) !- DO i=1,n CALL random_number (rndnum) ipique = INT(rndnum*n_rest)+1 ordre(i) = restant(ipique) restant(ipique:n_rest-1) = restant(ipique+1:n_rest) n_rest = n_rest-1 ENDDO !--------------------- END SUBROUTINE permute !- !=== !----------------- END MODULE weather