! MODULE readdim2 !--------------------------------------------------------------------- !- $Header: /home/ssipsl/CVSREP/sechiba/readdim2.f90,v 1.16 2000/05/11 14:48:27 ssipsl Exp $ !- USE ioipsl USE weather !- IMPLICIT NONE !- PRIVATE PUBLIC :: forcing_read, forcing_info, forcing_grid !- INTEGER, SAVE :: iim_full, jjm_full, llm_full, ttm_full INTEGER, SAVE :: iim_zoom, jjm_zoom, ist_zoom, jst_zoom REAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: data_full, lon_full, lat_full REAL, SAVE, ALLOCATABLE, DIMENSION(:) :: lev_full INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: itau, i_index, j_index INTEGER, SAVE :: i_test, j_test LOGICAL, SAVE :: allow_weathergen, weathergen, interpol REAL, SAVE :: limit_west, limit_east, limit_north, limit_south REAL, SAVE :: merid_res, zonal_res !- CONTAINS !- !===================================================================== !- SUBROUTINE forcing_info(filename, iim, jjm, llm, tm, date0, dt_force,& & force_id) !--------------------------------------------------------------------- ! !- This subroutine will get all the info from the forcing file and !- prepare for the zoom if needed. !- It returns to the caller the sizes of the data it will receive at !- the forcing_read call. This is important so that the caller can !- allocate the right space. !- !- filename : name of the file to be opened !- iim : size in x of the forcing data !- jjm : size in y of the forcing data !- llm : number of levels in the forcing data (not yet used) !- tm : Time dimension of the forcing !- date0 : The date at which the forcing file starts (julian days) !- dt-force : time-step of the forcing file in seconds !- force_id : ID of the forcing file !- !- ARGUMENTS !- CHARACTER(LEN=*) :: filename INTEGER :: iim, jjm, llm, tm REAL :: date0, dt_force INTEGER :: force_id !- LOCAL REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat REAL :: juld_1, juld_2 LOGICAL :: debug = .FALSE. !- CALL flininfo(filename, iim_full, jjm_full, llm_full, ttm_full, force_id) !- if ( debug ) WRITE(*,*) 'Details from forcing file :', iim_full, jjm_full, llm_full, ttm_full !- ALLOCATE(lon(iim_full, jjm_full), lat(iim_full, jjm_full)) ALLOCATE(itau(ttm_full)) ALLOCATE(data_full(iim_full, jjm_full),lon_full(iim_full, jjm_full),& & lat_full(iim_full, jjm_full)) ALLOCATE(lev_full(llm_full)) ALLOCATE(i_index(iim_full), j_index(jjm_full)) !- CALL flinopen & & (filename, .FALSE., iim_full, jjm_full, llm_full, lon_full, lat_full, & & lev_full, ttm_full, itau, date0, dt_force, force_id) !- !- What are the alowed options for the temportal interpolation !- !Config Key = ALLOW_WEATHERGEN !Config Desc = Allow weather generator to create data !Config Def = n !Config Help = This flag allows the forcing-reader to generate !Config synthetic data if the data in the file is too sparse !Config and the temporal resolution would not be enough to !Config run the model. !- allow_weathergen = .FALSE. CALL getin('ALLOW_WEATHERGEN',allow_weathergen) !- IF (ttm_full .GE. 2) THEN juld_1 = itau2date(itau(1), date0, dt_force) juld_2 = itau2date(itau(2), date0, dt_force) ELSE WRITE(*,*) 'forcing_info : What is that only one time step in the forcing file ' WRITE(*,*) 'That can not be right ' STOP ENDIF !- What is the distance between the two first states. From this we will deduce what is !- to be done. weathergen = .FALSE. interpol = .FALSE. !- IF ( ABS(ABS(juld_2-juld_1)-30.) .LE. 2.) THEN IF ( allow_weathergen ) THEN weathergen = .TRUE. WRITE(*,*) 'Using weather generator.' ELSE WRITE(*,*) 'This seems to be a monthly file.' WRITE(*,*) 'We should use a weather generator with this file.' WRITE(*,*) 'This should be allowed in the run.def' STOP ENDIF ELSEIF ( ABS(juld_1-juld_2) .LE. 1./4.) THEN interpol = .TRUE. WRITE(*,*) 'We will interpolate between the forcing data time steps.' ELSE ! Using the weather generator with data other than monthly ones probably ! needs some thinking. WRITE(*,*) 'We cannot do anything with these forcing data.' WRITE(*,*) 'The time step is not suitable:',ABS(juld_1-juld_2),' days.' STOP ENDIF !- !- redefine the forcing time step if the weather generator is activated !- IF ( weathergen ) THEN !Config Key = DT_WEATHGEN !Config Desc = Calling frequency of weather generator (s) !Config If = ALLOW_WEATHERGEN !Config Def = 1800. !Config Help = Determines how often the weather generator !Config is called (time step in s). Should be equal !Config to or larger than Sechiba's time step (say, !Config up to 6 times Sechiba's time step or so). dt_force = 1800. CALL getin('DT_WEATHGEN',dt_force) ENDIF !- !- Define the zoom !- !Config Key = LIMIT_WEST !Config Desc = Western limit of region !Config Def = -180. !Config Help = Western limit of the region we are !Config interested in. Between -180 and +180 degrees !Config The model will use the smalest regions from !Config region specified here and the one of the forcing file. !- limit_west = -180. CALL getin('LIMIT_WEST',limit_west) !- !Config Key = LIMIT_EAST !Config Desc = Eastern limit of region !Config Def = 180. !Config Help = Eastern limit of the region we are !Config interested in. Between -180 and +180 degrees !Config The model will use the smalest regions from !Config region specified here and the one of the forcing file. !- limit_east = 180. CALL getin('LIMIT_EAST',limit_east) !- !Config Key = LIMIT_NORTH !Config Desc = Northern limit of region !Config Def = 90. !Config Help = Northern limit of the region we are !Config interested in. Between +90 and -90 degrees !Config The model will use the smalest regions from !Config region specified here and the one of the forcing file. !- limit_north = 90. CALL getin('LIMIT_NORTH',limit_north) !- !Config Key = LIMIT_SOUTH !Config Desc = Southern limit of region !Config Def = -90. !Config Help = Southern limit of the region we are !Config interested in. Between 90 and -90 degrees !Config The model will use the smalest regions from !Config region specified here and the one of the forcing file. !- limit_south = -90. CALL getin('LIMIT_SOUTH',limit_south) !- !- Calculate domain size !- IF ( interpol ) THEN !- !- If we use temporal interpolation, then we cannot change the resolution (yet?) !- CALL domain_size (limit_west, limit_east, limit_north, limit_south,& & iim_full, jjm_full, lon_full, lat_full, iim_zoom, jjm_zoom,& & i_index, j_index) !- !- If we use the weather generator, then we read zonal and meridional resolutions !- This should be unified one day... !- ELSEIF ( weathergen ) THEN !- !Config Key = MERID_RES !Config Desc = North-South Resolution !Config Def = 2. !Config If = ALLOW_WEATHERGEN !Config Help = North-South Resolution of the region we are !Config interested in. In degrees !- merid_res = 2. CALL getin('MERID_RES',merid_res) !- !Config Key = ZONAL_RES !Config Desc = East-West Resolution !Config Def = 2. !Config If = ALLOW_WEATHERGEN !Config Help = East-West Resolution of the region we are !Config interested in. In degrees !- zonal_res = 2. CALL getin('ZONAL_RES',zonal_res) !- CALL weathgen_domain_size (limit_west,limit_east,limit_north,limit_south, & zonal_res,merid_res,iim_zoom,jjm_zoom) !- !- Number of time steps is meaningless in this case !- ttm_full = HUGE( ttm_full ) !- ELSE !- STOP 'ERROR: Neither interpolation nor weather generator is specified.' !- ENDIF !- !- Transfer the right information to the caller !- iim = iim_zoom jjm = jjm_zoom llm = 1 tm = ttm_full !- END SUBROUTINE forcing_info !- !- !===================================================================== SUBROUTINE forcing_read & & (filename, rest_id, lrstread, lrstwrite, & & itauin, istp, itau_split, split, nb_spread, netrad_cons, date0, & & dt_force, iim, jjm, lon, lat, ttm, & & swdown, precip, snowf, tair, & & u, v, qair, pb, lwdown, & & kindex, nbindex, force_id) !--------------------------------------------------------------------- !- filename : name of the file to be opened !- rest_id : ID of restart file !- lrstread : read restart file? !- lrstwrite : write restart file? !- itauin : time step for which we need the data !- istp : time step for restart file !- itau_split : Where are we within the splitting !- of the time-steps of the forcing files !- (it decides IF we READ or not) !- split : The number of time steps we do !- between two time-steps of the forcing !- nb_spread : Over how many time steps do we spread the precipitation !- netrad_cons: flag that decides if net radiation should be conserved. !- date0 : The date at which the forcing file starts (julian days) !- dt_force : time-step of the forcing file in seconds !- iim : Size of the grid in x !- jjm : size of the grid in y !- lon : Longitudes !- lat : Latitudes !- ttm : number of time steps in all in the forcing file !- swdown : Downward solar radiation (W/m^2) !- precip : Precipitation (Rainfall) (kg/m^2s) !- snowf : Snowfall (kg/m^2s) !- tair : 2m air temperature (K) !- u and v : 2m (in theory !) wind speed (m/s) !- qair : 2m humidity (kg/kg) !- pb : Surface pressure (Pa) !- lwdown : Downward long wave radiation (W/m^2) !- kindex : Index of all land-points in the data !- (used for the gathering) !- nbindex : Number of land points !- force_id : FLINCOM file id. !- It is used to close the file at the end of the run. !--------------------------------------------------------------------- IMPLICIT NONE !- CHARACTER(LEN=*) :: filename INTEGER :: rest_id LOGICAL :: lrstread, lrstwrite INTEGER :: itauin, istp, itau_split, split, nb_spread LOGICAL :: netrad_cons REAL :: date0, dt_force INTEGER :: iim, jjm, ttm REAL,DIMENSION(iim,jjm) :: & & lon, lat, swdown, precip, snowf, tair, u, v, qair, pb, lwdown INTEGER,DIMENSION(iim*jjm) :: kindex INTEGER :: nbindex, force_id !- INTEGER :: i,j REAL, SAVE :: un_an, un_jour = 0.0 !- !- Initialize if necessary ! IF ( un_jour .LE. 0 ) THEN CALL ioget_calendar(un_an, un_jour) ENDIF ! IF ( interpol ) THEN !- CALL forcing_read_interpol & (filename, itauin, itau_split, split, nb_spread, netrad_cons, date0, & dt_force, iim, jjm, lon, lat, ttm, swdown, precip, snowf, tair, & u, v, qair, pb, lwdown, kindex, nbindex, force_id) !- ELSEIF ( weathergen ) THEN !- CALL weathgen_main (itauin, istp, filename, force_id, iim, jjm, 1, & rest_id, lrstread, lrstwrite, & limit_west, limit_east, limit_north, limit_south, & zonal_res, merid_res, lon, lat, date0, & kindex, nbindex, dt_force, & swdown, precip, snowf, tair, u, v, qair, pb, lwdown) !- !- Transform from our old convention to the ALMA one. !- (Temporary fix) !- DO i=1,iim DO j=1,jjm precip(i,j) = precip(i,j)/un_jour snowf(i,j) = snowf(i,j)/un_jour tair(i,j) = tair(i,j)-273.15 pb(i,j) = pb(i,j)*100. ENDDO ENDDO !- ELSE !- STOP 'ERROR: Neither interpolation nor weather generator is specified.' !- ENDIF ! !------------------------- END SUBROUTINE forcing_read !===================================================================== !- !- !===================================================================== SUBROUTINE forcing_read_interpol & & (filename, itauin, itau_split, split, nb_spread, netrad_cons, date0, & & dt_force, iim, jjm, lon, lat, ttm, swdown, precip, snowf, tair, & & u, v, qair, pb, lwdown, kindex, nbindex, force_id) !--------------------------------------------------------------------- !- filename : name of the file to be opened !- itauin : time step for which we need the data !- itau_split : Where are we within the splitting !- of the time-steps of the forcing files !- (it decides IF we READ or not) !- split : The number of time steps we do !- between two time-steps of the forcing !- nb_spread : Over how many time steps do we spread the precipitation !- netrad_cons: flag that decides if net radiation should be conserved. !- date0 : The date at which the forcing file starts (julian days) !- dt_force : time-step of the forcing file in seconds !- iim : Size of the grid in x !- jjm : size of the grid in y !- lon : Longitudes !- lat : Latitudes !- ttm : number of time steps in all in the forcing file !- swdown : Downward solar radiation (W/m^2) !- precip : Precipitation (kg/m^2s) !- snowf : Snow fall (kg/m^2s) !- tair : 2m air temperature (K) !- u and v : 2m (in theory !) wind speed (m/s) !- qair : 2m humidity (kg/kg) !- pb : Surface pressure (Pa) !- lwdown : Downward long wave radiation (W/m^2) !- kindex : Index of all land-points in the data !- (used for the gathering) !- nbindex : Number of land points !- force_id : FLINCOM file id. !- It is used to close the file at the end of the run. !--------------------------------------------------------------------- IMPLICIT NONE !- INTEGER,PARAMETER :: lm=1 !- CHARACTER(LEN=*) :: filename INTEGER :: itauin, itau_split, split, nb_spread LOGICAL :: netrad_cons REAL :: date0, dt_force INTEGER :: iim, jjm, ttm REAL,DIMENSION(iim,jjm) :: & & lon, lat, swdown, precip, snowf, tair, u, v, qair, pb, lwdown INTEGER,DIMENSION(iim*jjm) :: kindex INTEGER :: nbindex, force_id !- !! INTEGER :: tiim, tjjm, tllm INTEGER, SAVE :: last_read=0, last_split INTEGER, SAVE :: itau_read, itau_read_n=0, itau_read_np1=0 REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: & & swdown_n, precip_n, snowf_n, tair_n, u_n, v_n, qair_n, pb_n, lwdown_n, & & swdown_np1, precip_np1, snowf_np1, tair_np1, u_np1, v_np1, qair_np1, & & pb_np1, lwdown_np1, mean_sinang, sinang REAL, SAVE :: julian_for, un_jour, un_an REAL :: julian, jur, ss, rw INTEGER :: yy, mm, dd, is LOGICAL, SAVE :: initialized = .FALSE. LOGICAL :: check=.FALSE. !--------------------------------------------------------------------- IF (check) THEN WRITE(*,*) & & " ISLSCP READ : itau_read, itau_split : ",itauin,itau_split ENDIF !- itau_read = itauin !- !- This part initializes the reading of the forcing. As you can see !- we only go through here if both time steps are zero. !- IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN !- IF (check) WRITE(*,*) 'ALLOCATE all the memory needed' !- ALLOCATE & & (swdown_n(iim,jjm), precip_n(iim,jjm), snowf_n(iim,jjm), & & tair_n(iim,jjm), u_n(iim,jjm), v_n(iim,jjm), qair_n(iim,jjm), & & pb_n(iim,jjm), lwdown_n(iim,jjm)) ALLOCATE & & (swdown_np1(iim,jjm), precip_np1(iim,jjm), snowf_np1(iim,jjm), & & tair_np1(iim,jjm), u_np1(iim,jjm), v_np1(iim,jjm), qair_np1(iim,jjm), & & pb_np1(iim,jjm), lwdown_np1(iim,jjm)) ALLOCATE & & (mean_sinang(iim,jjm), sinang(iim,jjm)) !- IF (check) WRITE(*,*) 'Memory ALLOCATED' !- !- We need for the driver surface air temperature and humidity before the !- the loop starts. Thus we read it here. !- CALL flinget (force_id,'Tair',iim_full, jjm_full, 1, ttm, 1, 1, data_full) CALL forcing_zoom(data_full, tair) CALL flinget (force_id,'PSurf',iim_full, jjm_full, 0, ttm, 1, 1, data_full) CALL forcing_zoom(data_full, pb) CALL flinget (force_id,'Qair',iim_full, jjm_full, 1, ttm, 1, 1, data_full) CALL forcing_zoom(data_full, qair) !--- IF (check) WRITE(*,*) 'tair, psol and qair have been extracted' !--- !--- Create the index table !--- CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test) CALL ioget_calendar(un_an, un_jour) !--- ENDIF !--- IF (check) THEN WRITE(*,*) & & 'The dates : ',itau_read,itau_split,itau_read_n,itau_read_np1 ENDIF !--- !--- Here we do the work in case only interpolation is needed. !--- IF ( initialized .AND. interpol ) THEN !--- IF (itau_read+1 /= last_read) THEN !--- IF (itau_read_n == 0) THEN CALL flinget (force_id,'SWdown', iim_full, jjm_full, 0,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, swdown_np1) CALL flinget (force_id,'Rainf', iim_full, jjm_full, 0,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, precip_np1) CALL flinget (force_id,'Snowf', iim_full, jjm_full, 0,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, snowf_np1) CALL flinget (force_id,'Tair', iim_full, jjm_full, 1,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, tair_np1) CALL flinget (force_id,'Wind_N', iim_full, jjm_full, 1,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, u_np1) CALL flinget (force_id,'Wind_E', iim_full, jjm_full, 1,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, v_np1) CALL flinget (force_id,'Qair', iim_full, jjm_full, 1,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, qair_np1) CALL flinget (force_id,'PSurf', iim_full, jjm_full, 0,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, pb_np1) CALL flinget (force_id,'LWdown', iim_full, jjm_full, 0,ttm, & & itau_read,itau_read, data_full) CALL forcing_zoom(data_full, lwdown_np1) itau_read_n = itau_read itau_read_np1 = itau_read ENDIF !----- swdown_n(:,:) = swdown_np1(:,:) precip_n(:,:) = precip_np1(:,:) tair_n(:,:) = tair_np1(:,:) u_n(:,:) = u_np1(:,:) v_n(:,:) = v_np1(:,:) qair_n(:,:) = qair_np1(:,:) pb_n(:,:) = pb_np1(:,:) lwdown_n(:,:) = lwdown_np1(:,:) itau_read_n = itau_read_np1 !----- itau_read_np1 = itau_read+1 IF (itau_read_np1 > ttm) THEN WRITE(*,*) 'WARNING --WARNING --WARNING --WARNING ' WRITE(*,*) & & 'WARNING : We are going back to the start of the file' itau_read_np1 =1 ENDIF !----- !----- Get a reduced julian day ! !----- This is needed because we lack the precision on 32 bit machines. !----- julian_for = itau2date(itau_read-1, date0, dt_force) CALL ju2ymds (julian_for, yy, mm, dd, ss) CALL ymds2ju (yy,1,1,0.0, jur) julian_for = julian_for-jur !----- IF (check) THEN WRITE(*,*) 'Forcing for Julian day ',julian_for,'is read' ENDIF !----- CALL flinget (force_id,'SWdown', iim_full, jjm_full, 0,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, swdown_np1) CALL flinget (force_id,'Rainf', iim_full, jjm_full, 0,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, precip_np1) CALL flinget (force_id,'Snowf', iim_full, jjm_full, 0,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, snowf_np1) CALL flinget (force_id,'Tair', iim_full, jjm_full, 1,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, tair_np1) CALL flinget (force_id,'Wind_N', iim_full, jjm_full, 1,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, u_np1) CALL flinget (force_id,'Wind_E', iim_full, jjm_full, 1,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, v_np1) CALL flinget (force_id,'Qair', iim_full, jjm_full, 1,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, qair_np1) CALL flinget (force_id,'PSurf', iim_full, jjm_full, 0,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, pb_np1) CALL flinget (force_id,'LWdown', iim_full, jjm_full, 0,ttm, & & itau_read_np1,itau_read_np1, data_full) CALL forcing_zoom(data_full, lwdown_np1) !----- last_read = itau_read_np1 !----- !----- Compute mean solar angle for the comming period !----- IF (check) WRITE(*,*) 'Going into solarang', split, un_jour !----- mean_sinang(:,:) = 0.0 DO is=1,split julian = julian_for+(FLOAT(is-1)/split)*dt_force/un_jour CALL solarang (julian, iim, jjm, lon, lat, sinang) mean_sinang(:,:) = mean_sinang(:,:)+sinang(:,:) ENDDO mean_sinang(:,:) = mean_sinang(:,:)/split !----- ENDIF !--- !--- Do the interpolation IF (check) WRITE(*,*) 'Doing the interpolation between time steps' !--- rw = REAL(itau_split-1) qair(:,:) = ((qair_np1(:,:)-qair_n(:,:))/split)*rw+qair_n(:,:) tair(:,:) = ((tair_np1(:,:)-tair_n(:,:))/split)*rw+tair_n(:,:) pb(:,:) = ((pb_np1(:,:)-pb_n(:,:))/split)*rw+pb_n(:,:) u(:,:) = ((u_np1(:,:)-u_n(:,:))/split)*rw+u_n(:,:) v(:,:) = ((v_np1(:,:)-v_n(:,:))/split)*rw+v_n(:,:) !--- !--- Here we need to allow for an option !--- where radiative energy is conserved !--- IF ( netrad_cons ) THEN lwdown(:,:) = lwdown_n(:,:) ELSE lwdown(:,:) = ((lwdown_np1(:,:)-lwdown_n(:,:))/split)*rw & & +lwdown_n(:,:) ENDIF !--- !--- For the solar radiation we decompose the mean value !--- using the zenith angle of the sun if the time step in the forcing data is !---- more than an hour. Else we use the standard linera interpolation !---- IF (check) WRITE(*,*) 'Ready to deal with the solar radiation' !---- IF ( dt_force .GT. 3600. ) THEN !--- IF ( netrad_cons ) THEN WRITE(*,*) 'Solar radiation can not be conserved with a timestep of ', dt_force ENDIF !--- julian = julian_for+(rw/split)*dt_force/un_jour IF (check) THEN WRITE(*,'(a,f20.10,2I3)') & & 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split ENDIF !--- CALL solarang(julian, iim, jjm, lon, lat, sinang) !--- WHERE (mean_sinang > 0.) swdown = swdown_n *sinang/mean_sinang ELSEWHERE swdown = 0.0 END WHERE !--- WHERE (swdown > 2000. ) swdown = 2000. END WHERE ELSE !--- IF ( netrad_cons ) THEN swdown(:,:) = swdown_n(:,:) ELSE swdown(:,:) = ((swdown_np1(:,:)-swdown_n(:,:))/split)*rw & & +swdown_n(:,:) ENDIF !--- ENDIF IF (check) THEN WRITE(*,*) 'SWNET at ',itau_split,' : ',swdown_n(i_test, j_test), & & ' < ', swdown(i_test, j_test), ' < ', swdown_np1(i_test, j_test) WRITE(*,*) 'TAIR at ',itau_split,' :',tair_n(i_test, j_test), & & ' < ', tair(i_test, j_test), ' < ', tair_np1(i_test, j_test) WRITE(*,*) 'QAIR at ',itau_split,' :',qair_n(i_test, j_test), & & ' < ', qair(i_test, j_test), ' < ', qair_np1(i_test, j_test) ENDIF !--- !--- For precip we suppose that the rain !--- is the summ over the next 6 hours !--- IF (itau_split <= nb_spread) THEN precip(:,:) = precip_n(:,:)/nb_spread ELSE precip(:,:) = 0.0 ENDIF !--- ENDIF !--- !--- Here we might put the call to the weather generator ... one day. !--- Pour le moment, le branchement entre interpolation et generateur de temps !--- est fait au-dessus. !--- !- IF ( initialized .AND. weathergen ) THEN !- .... !- ENDIF !--- !--- At this point the code should be initialized. If not we have a problem ! !--- IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN !--- initialized = .TRUE. !--- ELSE IF ( .NOT. initialized ) THEN WRITE(*,*) 'Why is the code forcing_read not initialized ?' WRITE(*,*) 'Have you called it with both time-steps set to zero ?' STOP ENDIF ENDIF ! !------------------------- END SUBROUTINE forcing_read_interpol !===================================================================== !- SUBROUTINE forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test) !--- !--- This subroutine finds the indices of the land points over which the land !--- surface scheme is going to run. !--- !- ARGUMENTS INTEGER :: iim, jjm, nbindex, i_test, j_test LOGICAL :: check INTEGER :: kindex(nbindex) REAL :: tair(iim,jjm) !- !- LOCAL INTEGER :: i, j, ig !- !- ig = 0 i_test = 0 j_test = 0 !--- IF (MINVAL(tair(:,:)) < 100.) THEN !----- In this case the 2m temperature is in Celsius DO j=1,jjm DO i=1,iim IF (tair(i,j) < 100.) THEN ig = ig+1 kindex(ig) = (j-1)*iim+i ! ! Here we find at random a land-point on which we can do ! some printouts for test. ! IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN i_test = i j_test = j IF (check) THEN WRITE(*,*) 'The test point chosen for output is : ', i_test, j_test ENDIF ENDIF ENDIF ENDDO ENDDO ELSE !----- 2m temerature is in Kelvin DO j=1,jjm DO i=1,iim IF (tair(i,j) < 500.) THEN ig = ig+1 kindex(ig) = (j-1)*iim+i ENDIF ENDDO ENDDO ENDIF !--- nbindex = ig !--- END SUBROUTINE forcing_landind !- !===================================================================== !- SUBROUTINE forcing_grid(iim,jjm,llm,lon,lat,lev) !- !- This subroutine calculates the longitudes and latitudes of the model grid. !- INTEGER, INTENT(in) :: iim, jjm, llm REAL, DIMENSION(iim,jjm), INTENT(out) :: lon, lat REAL, DIMENSION(llm), INTENT(out) :: lev !- INTEGER :: i,j REAL :: zlev !- !- Should be unified one day !- IF ( weathergen ) THEN !- DO i = 1, iim lon(i,:) = limit_west + zonal_res/2. + & FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim) ENDDO !- DO j = 1, jjm lat(:,j) = limit_north - merid_res/2. - & FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm) ENDDO !- !Config Key = HEIGHT_LEV1 !Config If = ALLOW_WEATHERGEN !Config Desc = Height at which the atmospheric variable are given !Config Def = 2.0 !Config Help = The atmospheric variables (temperature and specific !Config humidity are measured at a specific level. !Config The height of this level is needed to compute !Config corectly the turbulent transfer coefficients. !Config Look at the description of the forcing !Config DATA for the correct value. !- zlev = 2.0 CALL getin('HEIGHT_LEV1', zlev) lev(:) = zlev !- ELSEIF ( interpol ) THEN !- CALL forcing_zoom(lon_full, lon) CALL forcing_zoom(lat_full, lat) lev(:) = lev_full(:) !- ELSE !- STOP 'Neither weather generator nor temporal interpolation is specified.' !- ENDIF !- END SUBROUTINE forcing_grid !- !===================================================================== !- SUBROUTINE forcing_zoom(x_f, x_z) !- !- This subroutine takes the slab of data over which we wish to run the model. !- REAL :: x_f(iim_full, jjm_full), x_z(iim_zoom, jjm_zoom) !- INTEGER :: i,j !- DO i=1,iim_zoom DO j=1,jjm_zoom x_z(i,j) = x_f(i_index(i),j_index(j)) ENDDO ENDDO !- END SUBROUTINE forcing_zoom !- !===================================================================== SUBROUTINE solarang (julian, iim, jjm, lon, lat, csang) !--------------------------------------------------------------------- !- This subroutine computes the solar angle according to the method !- used by GSWP and developed by J.C. Morill. !- See for further details : !- http://www.atmo.arizona.edu/~morrill/swrad.html !--------------------------------------------------------------------- USE calendar !- INTEGER,INTENT(in) :: iim, jjm REAL,INTENT(in) :: julian REAL,DIMENSION(iim,jjm), INTENT(in) :: lon, lat REAL,DIMENSION(iim,jjm), INTENT(out) :: csang !- REAL :: un_an, un_jour, pi, gamma, dec, decd REAL :: et, gmt, le, ls, lcorr, latime, omega, omegad REAL :: llatd, llat INTEGER :: igmt, ilon, ilat INTEGER,SAVE,ALLOCATABLE :: zone(:), lhour(:) ! LOGICAL :: check = .FALSE. !--------------------------------------------------------------------- ! pi = 4.*ATAN(1.) IF (check) WRITE(*,*) 'We get the right calendar information' CALL ioget_calendar(un_an, un_jour) !- !- 1) Day angle gamma !- gamma = 2.*pi*MOD(julian,un_an)/un_an !- !- 2) Solar declination (assumed constant for a 24 hour period) page 7 !- in radians !- IF (check) WRITE(*,*) 'Solar declination' ! dec = ( 0.006918-0.399912*COS(gamma)+0.070257*SIN(gamma) & & -0.006758*COS(2*gamma)+0.000907*SIN(2*gamma) & & -0.002697*COS(3*gamma)+0.00148*SIN(3*gamma)) decd = dec*(180/pi) !- !- maximum error 0.0006 rad (<3'), leads to error !- of less than 1/2 degree in zenith angle !- IF (check) WRITE(*,*) 'Equation of time' !- 3) Equation of time page 11 !- et = ( 0.000075+0.001868*COS(gamma)-0.032077*SIN(gamma)& & -0.014615*COS(2*gamma)-0.04089*SIN(2*gamma))*229.18 !- !- Get the time zones for the current time !- gmt = 24.*(julian-INT(julian)) IF (.NOT.ALLOCATED(zone)) ALLOCATE(zone(iim)) IF (.NOT.ALLOCATED(lhour)) ALLOCATE(lhour(iim)) !- igmt = NINT(gmt) IF (check) WRITE(*,*) 'Get time zone' CALL time_zone(igmt, iim, jjm, lon, zone, lhour) !- !- Start a loop through the grid !- IF (check) WRITE(*,*) 'Start a loop through the grid' DO ilon=1,iim !--- !--- 4) Local apparent time page 13 !--- !--- ls standard longitude (nearest 15 degree meridian) !--- le local longtitude !--- lhour local standard time !--- latime local apparent time !--- lcorr longitudunal correction (minutes) !--- le = lon(ilon,1) ls = ((zone(ilon)-1)*15)-180. lcorr = 4.*(ls-le)*(-1) latime = lhour(ilon)+lcorr/60.+et/60. IF (latime < 0.) latime = latime+24 IF (latime > 24.) latime = latime-24 !--- !--- 5) Hour angle omega page 15 !--- !--- hour angle is zero at noon, postive in the morning !--- It ranges from 180 to -180 !--- !--- omegad is omega in degrees, omega is in radians !--- omegad = (latime-12.)*(-15.) omega = omegad*pi/180. !--- DO ilat=1,jjm !----- !----- 6) Zenith angle page 15 !----- !----- csang cosine of zenith angle (radians) !----- llatd = local latitude in degrees !----- llat = local latitude in radians !----- llatd = lat(1,ilat) llat = llatd*pi/180. csang(ilon,ilat) = & & MAX(0.,SIN(dec)*SIN(llat)+COS(dec)*COS(llat)*COS(omega)) ENDDO ENDDO !---------------------- END SUBROUTINE solarang !===================================================================== SUBROUTINE time_zone (gmt, iim, jjm, lon, zone, lhour) !--------------------------------------------------------------------- INTEGER :: gmt, iim, jjm, zone(iim), lhour(iim) REAL :: lon(iim,jjm) !- INTEGER :: deg, i, ilon, change !--------------------------------------------------------------------- DO ilon=1,iim !--- !--- Convert longitude index to longtitude in degrees !--- deg = lon(ilon,1) !--- !--- Determine into which time zone (15 degree interval) the !--- longitude falls. !--- DO i=1,25 IF (deg < (-187.5+(15*i))) THEN zone(ilon) = i IF (zone(ilon) == 25) zone(ilon) = 1 EXIT ENDIF ENDDO !--- !--- Calculate change (in number of hours) from GMT time to !--- local hour. Change will be negative for zones < 13 and !--- positive for zones > 13. !--- !--- There is also a correction for lhour < 0 and lhour > 23 !--- to lhour between 0 and 23. !--- IF (zone(ilon) < 13) THEN change = zone(ilon)-13 lhour(ilon) = gmt+change ENDIF !--- IF (zone(ilon) == 13) THEN lhour(ilon) = gmt ENDIF !--- IF (zone(ilon) > 13) THEN change = zone(ilon)-13 lhour(ilon) = gmt+change ENDIF IF (lhour(ilon) < 0) lhour(ilon) = lhour(ilon)+24 IF (lhour(ilon) > 23) lhour(ilon) = lhour(ilon)-24 !--- ENDDO !----------------------- END SUBROUTINE time_zone !===================================================================== ! ! --------------------------------------------------------------------- ! SUBROUTINE domain_size (limit_west, limit_east, limit_north, limit_south, & & iim_f, jjm_f, lon, lat, iim, jjm, iind, jind) IMPLICIT NONE ! ! ARGUMENTS ! REAL, INTENT(inout) :: limit_west,limit_east,limit_north,limit_south INTEGER, INTENT(in) :: iim_f, jjm_f REAL, INTENT(in) :: lon(iim_f, jjm_f), lat(iim_f, jjm_f) INTEGER, INTENT(out) :: iim,jjm INTEGER, INTENT(out) :: iind(iim_f), jind(jjm_f) ! ! LOCAL ! INTEGER :: i, j REAL :: lolo LOGICAL :: over_dateline = .FALSE. ! ! IF ( ( ABS(limit_east) .GT. 180. ) .OR. & ( ABS(limit_west) .GT. 180. ) ) THEN WRITE(*,*) 'PROBLEME longitudes.' WRITE(*,*) 'Limites Ouest, Est: ',limit_west,limit_east STOP ENDIF ! IF ( limit_west .GT. limit_east ) over_dateline = .TRUE. ! IF ( ( limit_south .LT. -90. ) .OR. & ( limit_north .GT. 90. ) .OR. & ( limit_south .GE. limit_north ) ) THEN WRITE(*,*) 'PROBLEME latitudes.' WRITE(*,*) 'Limites Nord, Sud: ',limit_north,limit_south STOP ENDIF ! ! Here we assume that the grid of the forcing data is regular. Else we would have ! to do more work to find the index table. ! iim = 0 DO i=1,iim_f ! lolo = lon(i,1) IF ( lon(i,1) .GT. 180. ) lolo = lon(i,1) - 360. IF ( lon(i,1) .LT. -180. ) lolo = lon(i,1) + 360. ! IF ( over_dateline ) THEN IF ( lolo .LT. limit_west .OR. lolo .GT. limit_east ) THEN iim = iim + 1 iind(iim) = i ENDIF ELSE IF ( lolo .GT. limit_west .AND. lolo .LT. limit_east ) THEN iim = iim + 1 iind(iim) = i ENDIF ENDIF ! ENDDO ! jjm = 0 DO j=1,jjm_f IF ( lat(1,j) .GT. limit_south .AND. lat(1,j) .LT. limit_north) THEN jjm = jjm + 1 jind(jjm) = j ENDIF ENDDO ! WRITE(*,*) 'Domain size: iim, jjm = ', iim, jjm END SUBROUTINE domain_size !------------------ END MODULE readdim2