!! !! This subroutine is the interface between Fortran 77 main program !! and SECHIBA. !! - Input fields are gathered to keep just continental points !! - call sechiba_main That's SECHIBA process. !! - Output fields are scattered to complete global fields !! !! @call sechiba_main !! @Version : $Revision: 1.9 $, $Date: 2000/04/05 14:48:08 $ !! !! @author Marie-Alice Foujols and Jan Polcher !! !! $Header: /home/ssipsl/CVSREP/sechiba/sechiba/intersurf.f90,v 1.9 2000/04/05 14:48:08 ssipsl Exp $ !! !f90doc MODULE intersurf MODULE intersurf USE IOIPSL USE reqdprec USE sechiba USE constantes USE constantes_veg IMPLICIT NONE PRIVATE PUBLIC :: intersurf_main INTERFACE intersurf_main MODULE PROCEDURE intersurf_main_2d, intersurf_main_1d END INTERFACE ! ! Global variables ! INTEGER(long),PARAMETER :: max_hist_level = 10 ! LOGICAL, SAVE :: l_first_intersurf=.TRUE. !! Initialisation has to be done one time ! INTEGER(long), SAVE :: hist_id, rest_id !! IDs for history and restart files INTEGER(long), SAVE :: hist_id_stom, rest_id_stom !! Dito for STOMATE ! INTEGER(long), SAVE :: itau_offset !! This offset is used to phase the !! calendar of the GCM or the driver. ! TYPE(control_type), SAVE :: control_flags !! Flags that (de)activate parts of the model ! REAL(stnd), ALLOCATABLE, DIMENSION (:,:), SAVE :: lalo !! Geographical coordinates INTEGER(long), ALLOCATABLE, DIMENSION (:,:), SAVE :: neighbours !! neighoring grid points if land REAL(stnd), ALLOCATABLE, DIMENSION (:,:), SAVE :: resolution !! size in x an y of the grid (m) ! CONTAINS ! !f90doc CONTAINS ! SUBROUTINE intersurf_main_2d (kjit, iim, jjm, kjpindex, kindex, xrdt, & & lrestart_read, lrestart_write, coupling, lon, lat, date0, & ! First level conditions & zlev, zlflu, u, v, qair, temp_air, epot_air, ccanopy, & ! Variables for the implicit coupling & cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, & ! Rain, snow, radiation and surface pressure & precip_rain, precip_snow, lwdown, swnet, swdown, pb, & ! Output : Fluxes & vevapp, fluxsens, fluxlat, run_off_tot, & ! Surface temperatures and surface properties & tsol_rad, temp_sol_new, albedo, emis, z0) ! routines called : sechiba_main ! ! ! interface description for dummy arguments ! input scalar INTEGER(long),INTENT (in) :: kjit !! Time step number INTEGER(long),INTENT (in) :: iim, jjm !! Dimension of input fields INTEGER(long),INTENT (in) :: kjpindex !! Number of continental points REAL(stnd),INTENT (in) :: xrdt !! Time step in seconds LOGICAL, INTENT (in) :: lrestart_read !! Logical for _restart_ file to read LOGICAL, INTENT (in) :: lrestart_write!! Logical for _restart_ file to write' CHARACTER*10, INTENT (in) :: coupling !! Describes the type of coupling REAL(stnd), INTENT (in) :: date0 !! Date at which kjit = 0 ! input fields INTEGER(long),DIMENSION (kjpindex), INTENT (in) :: kindex !! Index for continental points REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: u !! Lowest level wind speed REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: v !! Lowest level wind speed REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: zlev !! Height of first layer REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: zlflu !! Height of between the two first fluxes REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: qair !! Lowest level specific humidity REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: precip_rain !! Rain precipitation REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: precip_snow !! Snow precipitation REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: lwdown !! Down-welling long-wave flux REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: swnet !! Net surface short-wave flux REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: swdown !! Downwelling surface short-wave flux REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: temp_air !! Air temperature in Kelvin REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: epot_air !! Air humidity REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: ccanopy !! CO2 concentration in the canopy REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: petAcoef !! Coeficients A from the PBL resolution REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: peqAcoef !! One for T and another for q REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: petBcoef !! Coeficients B from the PBL resolution REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: peqBcoef !! One for T and another for q REAL(stnd),DIMENSION (iim,jjm), INTENT(inout) :: cdrag !! Cdrag REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: pb !! Lowest level pressure REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: lon, lat !! Geographical coordinates ! output fields REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: z0 !! Surface roughness REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: run_off_tot !! Complete runoff REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: tsol_rad !! REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: vevapp !! Total of evaporation REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: temp_sol_new !! New soil temperature REAL(stnd),DIMENSION (iim,jjm,2), INTENT(out) :: albedo !! Albedo REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: fluxsens !! Sensible chaleur flux REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: fluxlat !! Latent chaleur flux REAL(stnd),DIMENSION (iim,jjm), INTENT(out) :: emis !! Emissivity ! LOCAL declaration ! work arrays to scatter and/or gather information just before/after sechiba_main call's ! and to keep output value for next call REAL(stnd),DIMENSION (kjpindex) :: zu !! Work array to keep u REAL(stnd),DIMENSION (kjpindex) :: zv !! Work array to keep v REAL(stnd),DIMENSION (kjpindex) :: zzlev !! Work array to keep zlev REAL(stnd),DIMENSION (kjpindex) :: zzlflu !! Work array to keep zlflu REAL(stnd),DIMENSION (kjpindex) :: zqair !! Work array to keep qair REAL(stnd),DIMENSION (kjpindex) :: zprecip_rain !! Work array to keep precip_rain REAL(stnd),DIMENSION (kjpindex) :: zprecip_snow !! Work array to keep precip_snow REAL(stnd),DIMENSION (kjpindex) :: zlwdown !! Work array to keep lwdown REAL(stnd),DIMENSION (kjpindex) :: zswnet !! Work array to keep swnet REAL(stnd),DIMENSION (kjpindex) :: zswdown !! Work array to keep swdown REAL(stnd),DIMENSION (kjpindex) :: ztemp_air !! Work array to keep temp_air REAL(stnd),DIMENSION (kjpindex) :: zepot_air !! Work array to keep epot_air REAL(stnd),DIMENSION (kjpindex) :: zccanopy !! Work array to keep ccanopy REAL(stnd),DIMENSION (kjpindex) :: zpetAcoef !! Work array to keep petAcoef REAL(stnd),DIMENSION (kjpindex) :: zpeqAcoef !! Work array to keep peqAcoef REAL(stnd),DIMENSION (kjpindex) :: zpetBcoef !! Work array to keep petBcoef REAL(stnd),DIMENSION (kjpindex) :: zpeqBcoef !! Work array to keep peqVcoef REAL(stnd),DIMENSION (kjpindex) :: zcdrag !! Work array to keep cdrag REAL(stnd),DIMENSION (kjpindex) :: zpb !! Work array to keep pb REAL(stnd),DIMENSION (kjpindex) :: zz0 !! Work array to keep z0 REAL(stnd),DIMENSION (kjpindex) :: zrun_off_tot !! Work array to keep run_off_tot REAL(stnd),DIMENSION (kjpindex) :: ztsol_rad !! Work array to keep tsol_rad REAL(stnd),DIMENSION (kjpindex) :: zvevapp !! Work array to keep vevapp REAL(stnd),DIMENSION (kjpindex) :: ztemp_sol_new !! Work array to keep temp_sol_new REAL(stnd),DIMENSION (kjpindex,2) :: zalbedo !! Work array to keep albedo REAL(stnd),DIMENSION (kjpindex) :: zfluxsens !! Work array to keep fluxsens REAL(stnd),DIMENSION (kjpindex) :: zfluxlat !! Work array to keep fluxlat REAL(stnd),DIMENSION (kjpindex) :: zemis !! Work array to keep emis INTEGER(long) :: i, j, ik INTEGER(long) :: itau_sechiba LOGICAL :: check = .FALSE. IF (l_first_intersurf) THEN ! IF ( check ) WRITE(*,*) 'Initialisation of intersurf' ! ! Configuration of SSL specific parameters ! CALL intsurf_config(control_flags) ! ! Create the internal coordinate table ! IF ( (.NOT.ALLOCATED(lalo))) THEN ALLOCATE(lalo(kjpindex,2)) ENDIF ! DO ik=1,kjpindex j = ((kindex(ik)-1)/iim) + 1 i = (kindex(ik) - (j-1)*iim) lalo(ik,1) = lat(i,j) lalo(ik,2) = lon(i,j) ENDDO !- !- Compute variable to help describe the grid !- once the points are gathered. !- IF ( (.NOT.ALLOCATED(neighbours))) THEN ALLOCATE(neighbours(kjpindex,4)) ENDIF IF ( (.NOT.ALLOCATED(resolution))) THEN ALLOCATE(resolution(kjpindex,2)) ENDIF ! CALL grid_stuff (kjpindex, iim, jjm, lon, lat, kindex, neighbours, resolution) ! CALL intsurf_restart(kjit, iim, jjm, lon, lat, date0, xrdt, control_flags, rest_id, rest_id_stom, itau_offset) ! CALL intsurf_history(iim, jjm, lon, lat, kjit+itau_offset, date0, xrdt, control_flags, hist_id, hist_id_stom) ! IF ( check ) WRITE(*,*) 'End of Initialisation of intersurf' ! ENDIF ! ! 1. gather input fields from kindex array ! Warning : I'm not sure this interface with one dimension array is the good one DO ik=1, kjpindex j = ((kindex(ik)-1)/iim) + 1 i = (kindex(ik) - (j-1)*iim) zu(ik) = u(i,j) zv(ik) = v(i,j) zzlev(ik) = zlev(i,j) zzlflu(ik) = zlflu(i,j) zqair(ik) = qair(i,j) zprecip_rain(ik) = precip_rain(i,j) zprecip_snow(ik) = precip_snow(i,j) zlwdown(ik) = lwdown(i,j) zswnet(ik) = swnet(i,j) zswdown(ik) = swdown(i,j) ztemp_air(ik) = temp_air(i,j) zepot_air(ik) = epot_air(i,j) zccanopy(ik) = ccanopy(i,j) zpetAcoef(ik) = petAcoef(i,j) zpeqAcoef(ik) = peqAcoef(i,j) zpetBcoef(ik) = petBcoef(i,j) zpeqBcoef(ik) = peqBcoef(i,j) zcdrag(ik) = cdrag(i,j) zpb(ik) = pb(i,j) ! ! keep output variables from previous sechiba calls (necessary for stack storage) ! zemis(ik) = emis(i,j) zalbedo(ik,1) = albedo(i,j,1) zalbedo(ik,2) = albedo(i,j,2) zz0(ik) = z0(i,j) ENDDO ! ! 3. call sechiba for continental points only ! IF ( check ) WRITE(*,*) 'Calling sechiba' ! ! Shift the time step to phase the two models ! itau_sechiba = kjit + itau_offset ! CALL sechiba_main (itau_sechiba, iim*jjm, kjpindex, kindex, xrdt, & & lrestart_read, lrestart_write, coupling, control_flags, & & lalo, neighbours, resolution, & ! First level conditions & zzlev, zzlflu, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, & ! Variables for the implicit coupling & zcdrag, zpetAcoef, zpeqAcoef, zpetBcoef, zpeqBcoef, & ! Rain, snow, radiation and surface pressure & zprecip_rain ,zprecip_snow, zlwdown, zswnet, zswdown, zpb, & ! Output : Fluxes & zvevapp, zfluxsens, zfluxlat, zrun_off_tot, & ! Surface temperatures and surface properties & ztsol_rad, ztemp_sol_new, zalbedo, zemis, zz0, & ! File ids & rest_id, hist_id, rest_id_stom, hist_id_stom ) ! IF ( check ) WRITE(*,*) 'out of SECHIBA' ! ! 4. scatter output fields ! DO ik=1, kjpindex j = ((kindex(ik)-1)/iim) + 1 i = (kindex(ik) - (j-1)*iim) z0(i,j) = zz0(ik) run_off_tot(i,j) = zrun_off_tot(ik) tsol_rad(i,j) = ztsol_rad(ik) vevapp(i,j) = zvevapp(ik) temp_sol_new(i,j) = ztemp_sol_new(ik) albedo(i,j,1) = zalbedo(ik,1) albedo(i,j,2) = zalbedo(ik,2) fluxsens(i,j) = zfluxsens(ik) fluxlat(i,j) = zfluxlat(ik) emis(i,j) = zemis(ik) cdrag(i,j) = zcdrag(ik) ENDDO ! IF ( .NOT. l_first_intersurf) THEN ! CALL histwrite (hist_id, 'evap', itau_sechiba, vevapp, iim*jjm, kindex) CALL histwrite (hist_id, 'temp_sol', itau_sechiba, temp_sol_NEW, iim*jjm, kindex) CALL histwrite (hist_id, 'tsol_max', itau_sechiba, temp_sol_NEW, iim*jjm, kindex) CALL histwrite (hist_id, 'tsol_min', itau_sechiba, temp_sol_NEW, iim*jjm, kindex) CALL histwrite (hist_id, 'fluxsens', itau_sechiba, fluxsens, iim*jjm, kindex) CALL histwrite (hist_id, 'swnet', itau_sechiba, swnet, iim*jjm, kindex) CALL histwrite (hist_id, 'alb_vis', itau_sechiba, albedo(:,:,1), iim*jjm, kindex) CALL histwrite (hist_id, 'alb_nir', itau_sechiba, albedo(:,:,2), iim*jjm, kindex) CALL histwrite (hist_id, 'tair', itau_sechiba, temp_air, iim*jjm, kindex) CALL histwrite (hist_id, 'qair', itau_sechiba, qair, iim*jjm, kindex) ! ENDIF ! l_first_intersurf = .FALSE. ! IF (long_print) WRITE (numout,*) ' intersurf_main done ' END SUBROUTINE intersurf_main_2d ! SUBROUTINE intersurf_main_1d (kjit, iim, jjm, kjpindex, kindex, xrdt, & & lrestart_read, lrestart_write, coupling, lon, lat, date0, & ! First level conditions & zlev, zlflu, u, v, qair, temp_air, epot_air, ccanopy, & ! Variables for the implicit coupling & cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, & ! Rain, snow, radiation and surface pressure & precip_rain, precip_snow, lwdown, swnet, swdown, pb, & ! Output : Fluxes & vevapp, fluxsens, fluxlat, run_off_tot, & ! Surface temperatures and surface properties & tsol_rad, temp_sol_new, albedo, emis, z0) ! routines called : sechiba_main ! ! ! interface description for dummy arguments ! input scalar INTEGER(long),INTENT (in) :: kjit !! Time step number INTEGER(long),INTENT (in) :: iim, jjm !! Dimension of input fields INTEGER(long),INTENT (in) :: kjpindex !! Number of continental points REAL(stnd),INTENT (in) :: xrdt !! Time step in seconds LOGICAL, INTENT (in) :: lrestart_read !! Logical for _restart_ file to read LOGICAL, INTENT (in) :: lrestart_write!! Logical for _restart_ file to write' CHARACTER*10, INTENT (in) :: coupling !! Describes the type of coupling REAL(stnd), INTENT (in) :: date0 !! Date at which kjit = 0 ! input fields INTEGER(long),DIMENSION (kjpindex), INTENT (in) :: kindex !! Index for continental points REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: u !! Lowest level wind speed REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: v !! Lowest level wind speed REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: zlev !! Height of first layer REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: zlflu !! Height of between the two first fluxes REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: qair !! Lowest level specific humidity REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: precip_rain !! Rain precipitation REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: precip_snow !! Snow precipitation REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: lwdown !! Down-welling long-wave flux REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: swnet !! Net surface short-wave flux REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: swdown !! Downwelling surface short-wave flux REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: temp_air !! Air temperature in Kelvin REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: epot_air !! Air humidity REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: ccanopy !! CO2 concentration in the canopy REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: petAcoef !! Coeficients A from the PBL resolution REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: peqAcoef !! One for T and another for q REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: petBcoef !! Coeficients B from the PBL resolution REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: peqBcoef !! One for T and another for q REAL(stnd),DIMENSION (iim*jjm), INTENT(inout) :: cdrag !! Cdrag REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: pb !! Lowest level pressure REAL(stnd),DIMENSION (iim*jjm), INTENT(in) :: lon, lat !! Geographical coordinates ! output fields REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: z0 !! Surface roughness REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: run_off_tot !! Complete runoff REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: tsol_rad !! REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: vevapp !! Total of evaporation REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: temp_sol_new !! New soil temperature REAL(stnd),DIMENSION (iim*jjm,2), INTENT(out) :: albedo !! Albedo REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: fluxsens !! Sensible chaleur flux REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: fluxlat !! Latent chaleur flux REAL(stnd),DIMENSION (iim*jjm), INTENT(out) :: emis !! Emissivity ! LOCAL declaration ! work arrays to scatter and/or gather information just before/after sechiba_main call's ! and to keep output value for next call REAL(stnd),DIMENSION (kjpindex) :: zu !! Work array to keep u REAL(stnd),DIMENSION (kjpindex) :: zv !! Work array to keep v REAL(stnd),DIMENSION (kjpindex) :: zzlev !! Work array to keep zlev REAL(stnd),DIMENSION (kjpindex) :: zzlflu !! Work array to keep zlflu REAL(stnd),DIMENSION (kjpindex) :: zqair !! Work array to keep qair REAL(stnd),DIMENSION (kjpindex) :: zprecip_rain !! Work array to keep precip_rain REAL(stnd),DIMENSION (kjpindex) :: zprecip_snow !! Work array to keep precip_snow REAL(stnd),DIMENSION (kjpindex) :: zlwdown !! Work array to keep lwdown REAL(stnd),DIMENSION (kjpindex) :: zswnet !! Work array to keep swnet REAL(stnd),DIMENSION (kjpindex) :: zswdown !! Work array to keep swdown REAL(stnd),DIMENSION (kjpindex) :: ztemp_air !! Work array to keep temp_air REAL(stnd),DIMENSION (kjpindex) :: zepot_air !! Work array to keep epot_air REAL(stnd),DIMENSION (kjpindex) :: zccanopy !! Work array to keep ccanopy REAL(stnd),DIMENSION (kjpindex) :: zpetAcoef !! Work array to keep petAcoef REAL(stnd),DIMENSION (kjpindex) :: zpeqAcoef !! Work array to keep peqAcoef REAL(stnd),DIMENSION (kjpindex) :: zpetBcoef !! Work array to keep petBcoef REAL(stnd),DIMENSION (kjpindex) :: zpeqBcoef !! Work array to keep peqVcoef REAL(stnd),DIMENSION (kjpindex) :: zcdrag !! Work array to keep cdrag REAL(stnd),DIMENSION (kjpindex) :: zpb !! Work array to keep pb REAL(stnd),DIMENSION (kjpindex) :: zz0 !! Work array to keep z0 REAL(stnd),DIMENSION (kjpindex) :: zrun_off_tot !! Work array to keep run_off_tot REAL(stnd),DIMENSION (kjpindex) :: ztsol_rad !! Work array to keep tsol_rad REAL(stnd),DIMENSION (kjpindex) :: zvevapp !! Work array to keep vevapp REAL(stnd),DIMENSION (kjpindex) :: ztemp_sol_new !! Work array to keep temp_sol_new REAL(stnd),DIMENSION (kjpindex,2) :: zalbedo !! Work array to keep albedo REAL(stnd),DIMENSION (kjpindex) :: zfluxsens !! Work array to keep fluxsens REAL(stnd),DIMENSION (kjpindex) :: zfluxlat !! Work array to keep fluxlat REAL(stnd),DIMENSION (kjpindex) :: zemis !! Work array to keep emis INTEGER(long) :: i, j, ik INTEGER(long) :: itau_sechiba REAL(stnd), ALLOCATABLE, DIMENSION(:,:) :: tmp_lon, tmp_lat LOGICAL :: check = .FALSE. IF (l_first_intersurf) THEN ! IF ( check ) WRITE(*,*) 'Initialisation of intersurf' ! ! Configuration of SSL specific parameters ! CALL intsurf_config(control_flags) ! ! Create the internal coordinate table ! IF ( (.NOT.ALLOCATED(tmp_lon))) THEN ALLOCATE(tmp_lon(iim,jjm)) ENDIF IF ( (.NOT.ALLOCATED(tmp_lat))) THEN ALLOCATE(tmp_lat(iim,jjm)) ENDIF ! IF ( (.NOT.ALLOCATED(lalo))) THEN ALLOCATE(lalo(kjpindex,2)) ENDIF ! DO i=1,iim DO j=1,jjm ik = (j-1)*iim + i tmp_lon(i,j) = lon(ik) tmp_lat(i,j) = lat(ik) ENDDO ENDDO ! lalo(:,1) = lat(:) lalo(:,2) = lon(:) !- !- Compute variable to help describe the grid !- once the points are gathered. !- IF ( (.NOT.ALLOCATED(neighbours))) THEN ALLOCATE(neighbours(kjpindex,4)) ENDIF IF ( (.NOT.ALLOCATED(resolution))) THEN ALLOCATE(resolution(kjpindex,2)) ENDIF ! CALL grid_stuff (kjpindex, iim, jjm, tmp_lon, tmp_lat, kindex, neighbours, resolution) ! CALL intsurf_restart(kjit, iim, jjm, tmp_lon, tmp_lat, date0, xrdt, control_flags, rest_id, rest_id_stom, itau_offset) ! CALL intsurf_history(iim, jjm, tmp_lon, tmp_lat, kjit+itau_offset, date0, xrdt, control_flags, hist_id, hist_id_stom) ! IF ( check ) WRITE(*,*) 'End of Initialisation of intersurf' ! ENDIF ! ! 1. gather input fields from kindex array ! Warning : I'm not sure this interface with one dimension array is the good one DO ik=1, kjpindex zu(ik) = u(kindex(ik)) zv(ik) = v(kindex(ik)) zzlev(ik) = zlev(kindex(ik)) zzlflu(ik) = zlflu(kindex(ik)) zqair(ik) = qair(kindex(ik)) zprecip_rain(ik) = precip_rain(kindex(ik)) zprecip_snow(ik) = precip_snow(kindex(ik)) zlwdown(ik) = lwdown(kindex(ik)) zswnet(ik) = swnet(kindex(ik)) zswdown(ik) = swdown(kindex(ik)) ztemp_air(ik) = temp_air(kindex(ik)) zepot_air(ik) = epot_air(kindex(ik)) zccanopy(ik) = ccanopy(kindex(ik)) zpetAcoef(ik) = petAcoef(kindex(ik)) zpeqAcoef(ik) = peqAcoef(kindex(ik)) zpetBcoef(ik) = petBcoef(kindex(ik)) zpeqBcoef(ik) = peqBcoef(kindex(ik)) zcdrag(ik) = cdrag(kindex(ik)) zpb(ik) = pb(kindex(ik)) ! ! keep output variables from previous sechiba calls (necessary for stack storage) ! zemis(ik) = emis(kindex(ik)) zalbedo(ik,1) = albedo(kindex(ik),1) zalbedo(ik,2) = albedo(kindex(ik),2) zz0(ik) = z0(kindex(ik)) ENDDO ! ! 3. call sechiba for continental points only ! IF ( check ) WRITE(*,*) 'Calling sechiba' ! ! Shift the time step to phase the two models ! itau_sechiba = kjit + itau_offset ! CALL sechiba_main (itau_sechiba, iim*jjm, kjpindex, kindex, xrdt, & & lrestart_read, lrestart_write, coupling, control_flags, & & lalo, neighbours, resolution, & ! First level conditions & zzlev, zzlflu, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, & ! Variables for the implicit coupling & zcdrag, zpetAcoef, zpeqAcoef, zpetBcoef, zpeqBcoef, & ! Rain, snow, radiation and surface pressure & zprecip_rain ,zprecip_snow, zlwdown, zswnet, zswdown, zpb, & ! Output : Fluxes & zvevapp, zfluxsens, zfluxlat, zrun_off_tot, & ! Surface temperatures and surface properties & ztsol_rad, ztemp_sol_new, zalbedo, zemis, zz0, & ! File ids & rest_id, hist_id, rest_id_stom, hist_id_stom ) ! IF ( check ) WRITE(*,*) 'out of SECHIBA' ! ! 4. scatter output fields ! DO ik=1, kjpindex z0(kindex(ik)) = zz0(ik) run_off_tot(kindex(ik)) = zrun_off_tot(ik) tsol_rad(kindex(ik)) = ztsol_rad(ik) vevapp(kindex(ik)) = zvevapp(ik) temp_sol_new(kindex(ik)) = ztemp_sol_new(ik) albedo(kindex(ik),1) = zalbedo(ik,1) albedo(kindex(ik),2) = zalbedo(ik,2) fluxsens(kindex(ik)) = zfluxsens(ik) fluxlat(kindex(ik)) = zfluxlat(ik) emis(kindex(ik)) = zemis(ik) cdrag(kindex(ik)) = zcdrag(ik) ENDDO ! IF ( .NOT. l_first_intersurf) THEN ! CALL histwrite (hist_id, 'evap', itau_sechiba, vevapp, iim*jjm, kindex) CALL histwrite (hist_id, 'temp_sol', itau_sechiba, temp_sol_NEW, iim*jjm, kindex) CALL histwrite (hist_id, 'tsol_max', itau_sechiba, temp_sol_NEW, iim*jjm, kindex) CALL histwrite (hist_id, 'tsol_min', itau_sechiba, temp_sol_NEW, iim*jjm, kindex) CALL histwrite (hist_id, 'fluxsens', itau_sechiba, fluxsens, iim*jjm, kindex) CALL histwrite (hist_id, 'swnet', itau_sechiba, swnet, iim*jjm, kindex) CALL histwrite (hist_id, 'alb_vis', itau_sechiba, albedo(:,1), iim*jjm, kindex) CALL histwrite (hist_id, 'alb_nir', itau_sechiba, albedo(:,2), iim*jjm, kindex) CALL histwrite (hist_id, 'tair', itau_sechiba, temp_air, iim*jjm, kindex) CALL histwrite (hist_id, 'qair', itau_sechiba, qair, iim*jjm, kindex) ! ENDIF ! l_first_intersurf = .FALSE. ! IF (long_print) WRITE (numout,*) ' intersurf_main done ' END SUBROUTINE intersurf_main_1d ! ! SUBROUTINE intsurf_config(control_flags) ! ! This subroutine reads all the configuration flags which control the behaviour of the model ! TYPE(control_type), INTENT(out) :: control_flags !! Flags that (de)activate parts of the model ! ! !Config Key = STOMATE_OK_CO2 !Config Desc = Activate CO2? !Config Def = n !Config Help = set to TRUE if photosynthesis is to be activated ! control_flags%ok_co2 = .FALSE. CALL getin('STOMATE_OK_CO2', control_flags%ok_co2) WRITE(*,*) 'photosynthesis: ', control_flags%ok_co2 ! !Config Key = STOMATE_OK_STOMATE !Config Desc = Activate STOMATE? !Config Def = n !Config Help = set to TRUE if STOMATE is to be activated ! control_flags%ok_stomate = .FALSE. CALL getin('STOMATE_OK_STOMATE',control_flags%ok_stomate) WRITE(*,*) 'STOMATE is activated: ',control_flags%ok_stomate ! !Config Key = STOMATE_OK_DGVM !Config Desc = Activate DGVM? !Config Def = n !Config Help = set to TRUE if DGVM is to be activated ! control_flags%ok_dgvm = .FALSE. CALL getin('STOMATE_OK_DGVM',control_flags%ok_dgvm) WRITE(*,*) 'LPJ is activated: ',control_flags%ok_dgvm ! !Config Key = STOMATE_WATCHOUT !Config Desc = STOMATE does minimum service !Config Def = n !Config Help = set to TRUE if you want STOMATE to read !Config and write its start files and keep track !Config of longer-term biometeorological variables. !Config This is useful if OK_STOMATE is not set, !Config but if you intend to activate STOMATE later. !Config In that case, this run can serve as a !Config spinup for longer-term biometeorological !Config variables. ! control_flags%stomate_watchout = .FALSE. CALL getin('STOMATE_WATCHOUT',control_flags%stomate_watchout) WRITE(*,*) 'STOMATE keeps an eye open: ',control_flags%stomate_watchout ! ! RETURN ! END SUBROUTINE intsurf_config ! ! ! SUBROUTINE intsurf_restart(istp, iim, jjm, lon, lat, date0, dt, control_flags, rest_id, rest_id_stom, itau_offset) ! ! This subroutine initialized the restart file for the land-surface scheme ! INTEGER(long), INTENT(in) :: istp !! Time step of the restart file INTEGER(long), INTENT(in) :: iim, jjm !! Size in x and y of the data to be handeled REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: lon, lat !! Logitude and latitude of the data points REAL(stnd) :: date0 !! The date at which itau = 0 REAL(stnd) :: dt !! Time step INTEGER(long), INTENT(out) :: rest_id, rest_id_stom !! ID of the restart file INTEGER(long), INTENT(out) :: itau_offset ! TYPE(control_type), INTENT(in) :: control_flags !! Flags that (de)activate parts of the model ! ! LOCAL ! CHARACTER(LEN=80) :: restname_in, restname_out, stom_restname_in, stom_restname_out REAL(stnd) :: dt_rest, date0_rest INTEGER(long) :: itau_dep INTEGER(long),PARAMETER :: llm=1 REAL(stnd), DIMENSION(llm) :: lev LOGICAL :: overwrite_time REAL(stnd) :: in_julian, rest_julian REAl(stnd) :: un_an, un_jour ! !Config Key = SECHIBA_restart_in !Config Desc = Name of restart to READ for initial conditions !Config Def = NONE !Config Help = This is the name of the file which will be opened !Config to extract the initial values of all prognostic !Config values of the model. This has to be a netCDF file. !Config Not truly COADS compliant. NONE will mean that !Config no restart file is to be expected. !- restname_in = 'NONE' CALL getin('SECHIBA_restart_in',restname_in) WRITE(*,*) 'INPUT RESTART_FILE', restname_in !- !Config Key = SECHIBA_rest_out !Config Desc = Name of restart files to be created by SECHIBA !Config Def = sechiba_rest_out.nc !Config Help = This variable give the name for !Config the restart files. The restart software within !Config IOIPSL will add .nc if needed. ! restname_out = 'restart_out.nc' CALL getin('SECHIBA_rest_out', restname_out) ! !Config Key = SECHIBA_reset_time !Config Desc = Option to overrides the time of the restart !Config Def = n !Config Help = This option allows the model to override the time !Config found in the restart file of SECHIBA with the time !Config of the first call. That is the restart time of the GCM. ! overwrite_time = .FALSE. CALL getin('SECHIBA_reset_time', overwrite_time) ! CALL ioget_calendar(un_an, un_jour) ! in_julian = itau2date(istp, date0, dt) date0_rest = date0 dt_rest = dt ! CALL restini( restname_in, iim, jjm, lon, lat, llm, lev, & & restname_out, itau_dep, date0_rest, dt_rest, rest_id, overwrite_time) ! ! itau_dep of SECHIBA is phased with the GCM if needed ! rest_julian = itau2date(itau_dep, date0_rest, dt_rest) ! IF ( ABS( in_julian - rest_julian) .GT. dt/un_jour ) THEN IF ( overwrite_time ) THEN WRITE(*,*) 'The SECHIBA restart is not for the same timestep as the GCM,' WRITE(*,*) 'the two are synchronized. The land-surface conditions can not impose' WRITE(*,*) 'the chronology of the simulation.' WRITE(*,*) 'Time step of the GCM :', istp, 'Julian day : ', in_julian WRITE(*,*) 'Time step of SECHIBA :', itau_dep, 'Julian day : ', rest_julian itau_offset = itau_dep - istp ELSE WRITE(*,*) 'IN -> OUT :', istp, '->', itau_dep WRITE(*,*) 'IN -> OUT :', in_julian, '->', rest_julian WRITE(*,*) 'SECHIBA''s restart file is not consistent with the one of the GCM' WRITE(*,*) 'Correct the time axis of the restart file or force the code to change it.' STOP ENDIF ELSE itau_offset = 0 ENDIF ! !===================================================================== !- 1.5 Restart file for STOMATE !===================================================================== IF ( control_flags%ok_stomate .OR. control_flags%stomate_watchout ) THEN !- ! STOMATE IS ACTIVATED !- !Config Key = STOMATE_RESTART_FILEIN !Config Desc = Name of restart to READ for initial conditions !Config of STOMATE !Config If = STOMATE_OK_STOMATE || STOMATE_WATCHOUT !Config Def = NONE !Config Help = This is the name of the file which will be opened !Config to extract the initial values of all prognostic !Config values of STOMATE. !- stom_restname_in = 'NONE' CALL getin('STOMATE_RESTART_FILEIN',stom_restname_in) WRITE(*,*) 'STOMATE INPUT RESTART_FILE', stom_restname_in !- !Config Key = STOMATE_RESTART_FILEOUT !Config Desc = Name of restart files to be created by STOMATE !Config If = STOMATE_OK_STOMATE || STOMATE_WATCHOUT !Config Def = stomate_restart.nc !Config Help = This is the name of the file which will be opened !Config to write the final values of all prognostic values !Config of STOMATE. !- stom_restname_out = 'stomate_restart.nc' CALL getin('STOMATE_RESTART_FILEOUT', stom_restname_out) WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE', stom_restname_out !- CALL restini (stom_restname_in, iim, jjm, lon, lat, llm, lev, & & stom_restname_out, itau_dep, date0_rest, dt_rest, rest_id_stom, overwrite_time) !- ENDIF ! END SUBROUTINE intsurf_restart SUBROUTINE intsurf_history(iim, jjm, lon, lat, istp_old, date0, dt, control_flags, hist_id, hist_id_stom) ! ! ! This subroutine initialized the history files for the land-surface scheme ! ! INTEGER(long), INTENT(in) :: iim, jjm !! Size in x and y of the data to be handeled REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: lon, lat !! Logitude and latitude of the data points INTEGER(long), INTENT(in) :: istp_old !! Time step counter REAL(stnd), INTENT(in) :: date0 !! Julian day at which istp=0 REAL(stnd), INTENT(in) :: dt !! Time step of the counter in seconds ! TYPE(control_type), INTENT(in) :: control_flags !! Flags that (de)activate parts of the model ! INTEGER(long), INTENT(out) :: hist_id, hist_id_stom !! History file identification for SECHIBA and STOMATE ! ! LOCAL ! CHARACTER(LEN=80) :: histname,stom_histname !! Name of history file REAL(stnd) :: dw !! frequency of history write (sec.) CHARACTER(LEN=30) :: flux_op !! Operations to be performed on fluxes CHARACTER(LEN=30) :: flux_sc !! Operations which do not include a scatter INTEGER(long) :: hist_level !! history output level (default is 10 => maximum output) CHARACTER(LEN=40),DIMENSION(max_hist_level) :: & & ave, avecels, avescatter, tmax, fluxop, fluxop_sc, & & tmincels, tmaxcels, once !! The various operation to be performed INTEGER(long) :: i INTEGER(long) :: hori_id !! ID of the default horizontal longitude and latitude map. INTEGER(long) :: vegax_id, solax_id, soltax_id !! ID's for two vertical coordinates INTEGER(long) :: hist_PFTaxis_id ! REAL(stnd),DIMENSION(nvm) :: veg REAL(stnd),DIMENSION(ngrnd) :: sol REAL(stnd),DIMENSION(nbsoil):: soltyp ! REAL(stnd) :: hist_days_stom !!- GK time step in days for this history file REAL(stnd) :: hist_dt_stom !!- GK time step in seconds for this history file INTEGER(long),PARAMETER :: npft = nvm-1 !!- GK Number of PFTs REAL(stnd),DIMENSION(npft) :: hist_PFTaxis !!- GK An axis we need for the history files ! !===================================================================== !- 3.0 Setting up the history files !===================================================================== !- 3.1 SECHIBA !===================================================================== !Config Key = OUTPUT_FILE !Config Desc = Name of file in which the output is going !Config to be written !Config Def = cabauw_out.nc !Config Help = This file is going to be created by the model !Config and will contain the output from the model. !Config This file is a truly COADS compliant netCDF file. !Config It will be generated by the hist software from !Config the IOIPSL package. !- histname='cabauw_out.nc' CALL getin('OUTPUT_FILE', histname) WRITE(*,*) 'OUTPUT_FILE', histname !- !Config Key = WRITE_STEP !Config Desc = Frequency in seconds at which to WRITE output !Config Def = 86400.0 !Config Help = This variables gives the frequency the output of !Config the model should be written into the netCDF file. !Config It does not affect the frequency at which the !Config operations such as averaging are done. !Config That is IF the coding of the calls to histdef !Config are correct ! !- dw = un_jour CALL getin('WRITE_STEP', dw) !- ! veg(1:nvm) = (/ (FLOAT(i),i=1,nvm) /) sol(1:ngrnd) = (/ (FLOAT(i),i=1,ngrnd) /) soltyp(1:nbsoil) = (/ (FLOAT(i),i=1,nbsoil) /) ! !- We need to flux averaging operation as when the data is written !- from within SECHIBA a scatter is needed. In the driver on the other !- hand the data is 2D and can be written is it is. !- WRITE(flux_op,'("ave(scatter(X*",F8.1,"))")') un_jour/dt WRITE(flux_sc,'("ave(X*",F8.1,")")') un_jour/dt WRITE(*,*) flux_op, un_jour/dt, dt, dw !- !Config Key = SECHIBA_HISTLEVEL !Config Desc = SECHIBA history output level (0..10) !Config Def = 5 !Config Help = Chooses the list of variables in the history file. !Config Values between 0: nothing is written; 10: everything is !Config written are available More details can be found on the web under documentation. !Config web under documentation. !- hist_level = 5 CALL getin('SECHIBA_HISTLEVEL', hist_level) !- WRITE(*,*) 'SECHIBA history level: ',hist_level IF ( (hist_level > max_hist_level).OR.(hist_level < 0) ) THEN STOP 'This history level is not allowed' ENDIF !- !- define operations as a function of history level. !- Above hist_level, operation='never' !- ave(1:10) = 'ave(X)' IF (hist_level < 10) THEN ave(hist_level+1:10) = 'never' ENDIF avecels(1:10) = 'ave(cels(X))' IF (hist_level < 10) THEN avecels(hist_level+1:10) = 'never' ENDIF avescatter(1:10) = 'ave(scatter(X))' IF (hist_level < 10) THEN avescatter(hist_level+1:10) = 'never' ENDIF tmincels(1:10) = 't_min(cels(X))' IF (hist_level < 10) THEN tmincels(hist_level+1:10) = 'never' ENDIF tmaxcels(1:10) = 't_max(cels(X))' IF (hist_level < 10) THEN tmaxcels(hist_level+1:10) = 'never' ENDIF tmax(1:10) = 't_max(X)' IF (hist_level < 10) THEN tmax(hist_level+1:10) = 'never' ENDIF fluxop(1:10) = flux_op IF (hist_level < 10) THEN fluxop(hist_level+1:10) = 'never' ENDIF fluxop_sc(1:10) = flux_sc IF (hist_level < 10) THEN fluxop_sc(hist_level+1:10) = 'never' ENDIF once(1:10) = 'once(scatter(X))' IF (hist_level < 10) THEN once(hist_level+1:10) = 'never' ENDIF !- CALL histbeg(histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & & istp_old, date0, dt, hori_id, hist_id) !- CALL histvert(hist_id, 'veget', 'Vegetation types', '-', & & nvm, veg, vegax_id) CALL histvert(hist_id, 'solth', 'Soil levels', 'm', & & ngrnd, sol, solax_id) CALL histvert(hist_id, 'soiltyp', 'Soil types', '-', & & nbsoil, soltyp, soltax_id) !- !- SECHIBA_HISTLEVEL = 1 !- CALL histdef(hist_id, 'evap', 'Evaporation', 'mm/d', & & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop_sc(1), dt,dw) CALL histdef(hist_id, 'temp_sol', 'Surface Temperature', 'C', & & iim,jjm, hori_id, 1,1,1, -99, 32, avecels(1), dt,dw) CALL histdef(hist_id, 'rain', 'Rainfall', 'mm/d', & & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(1), dt,dw) CALL histdef(hist_id, 'snowf', 'Snowfall', 'mm/d', & & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(1), dt,dw) CALL histdef(hist_id, 'netrad', 'Net radiation', 'W/m^2', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw) CALL histdef(hist_id, 'lai', 'Leaf Area Index', '-', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(1), dt,dw) !- !- SECHIBA_HISTLEVEL = 2 !- CALL histdef(hist_id, 'subli', 'Sublimation', 'mm/d', & & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(2), dt,dw) CALL histdef(hist_id, 'evapnu', 'Bare soil evaporation', 'mm/d', & & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(2), dt,dw) CALL histdef(hist_id, 'runoff', 'Surface runoff', 'mm/d', & & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(2), dt,dw) CALL histdef(hist_id, 'tair', 'Air Temperature', 'C', & & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw) CALL histdef(hist_id, 'qair', 'Air humidity', 'g/g', & & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw) CALL histdef(hist_id, 'alb_vis', 'Albedo visible', '-', & & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw) CALL histdef(hist_id, 'alb_nir', 'Albedo near infrared', '-', & & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw) CALL histdef(hist_id, 'z0', 'Surface roughness', '-', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(2), dt,dw) CALL histdef(hist_id, 'transpir', 'Transpiration', 'mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(2), dt,dw) CALL histdef(hist_id, 'inter', 'Interception loss', 'mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(2), dt,dw) !- !- SECHIBA_HISTLEVEL = 3 !- CALL histdef(hist_id, 'tsol_max', 'Maximum Surface Temperature',& & 'C', iim,jjm, hori_id, 1,1,1, -99, 32, tmaxcels(3), dt,dw) CALL histdef(hist_id, 'tsol_min', 'Minimum Surface Temperature',& & 'C', iim,jjm, hori_id, 1,1,1, -99, 32, tmincels(3), dt,dw) CALL histdef(hist_id, 'fluxsens', 'Sensible Heat Flux', 'W/m^2', & & iim,jjm, hori_id, 1,1,1, -99, 32, ave(3), dt,dw) CALL histdef(hist_id, 'snow', 'Snow mass', 'kg/m^2', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw) CALL histdef(hist_id, 'snowage', 'Snow age', '?', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw) CALL histdef(hist_id, 'snowice', 'Snow on ice', 'kg/m^2', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw) CALL histdef(hist_id, 'snowiceage', 'Snow age on ice', 'd', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw) CALL histdef(hist_id, 'vegetfrac', 'Fraction of vegetation', '-', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(3), dt,dw) !- !- SECHIBA_HISTLEVEL = 4 !- CALL histdef(hist_id, 'gqsb', 'Upper Soil Moisture', '-', & & iim,jjm, hori_id, 1, 1, 1, -99, 32, avescatter(4), dt,dw) CALL histdef(hist_id, 'bqsb', 'Lower Soil Moisture', '-', & & iim,jjm, hori_id, 1, 1, 1, -99, 32, avescatter(4), dt,dw) CALL histdef(hist_id, 'qsintveg', 'Water on canopy', 'Kg/m^2', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(4), dt,dw) CALL histdef(hist_id, 'rstruct', 'Structural resistance', 'mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(4), dt,dw) CALL histdef(hist_id, 'gpp', 'CO2 Assimilation', 'gC/m2/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(4), dt,dw) CALL histdef(hist_id, 'precisol', 'Throughfall', 'mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(4), dt,dw) CALL histdef(hist_id, 'hdry', 'Dry soil height', 'm', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(4), dt,dw) !- !- SECHIBA_HISTLEVEL = 5 !- CALL histdef(hist_id, 'swnet', 'Net solar radiation', 'W/m^2', & & iim,jjm, hori_id, 1,1,1, -99, 32, ave(5), dt,dw) CALL histdef(hist_id, 'temp_pheno', 'Temperature for Pheno', 'K', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(5), dt,dw) !- !- SECHIBA_HISTLEVEL = 6 !- CALL histdef(hist_id, 'ptn', 'Deep ground temperature', 'K', & & iim,jjm, hori_id, ngrnd, 1, ngrnd, solax_id, 32, ave(6), dt,dw) !- !- SECHIBA_HISTLEVEL = 7 !- !- !- SECHIBA_HISTLEVEL = 8 !- CALL histdef(hist_id, 'beta', 'Beta Function', '-', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw) CALL histdef(hist_id, 'qcdrag', 'Diffusion coeff.', '-', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw) CALL histdef(hist_id, 'vbeta1', 'Beta for sublimation', '-', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw) CALL histdef(hist_id, 'vbeta4', 'Beta for bare soil', '-', & & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw) CALL histdef(hist_id, 'vbetaco2', 'beta for CO2', 'mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(8), dt,dw) CALL histdef(hist_id, 'dtdt', 'Max value of dT/dt', 'K/s', & & iim,jjm, hori_id, 1,1,1, -99, 32, tmax(8), dt,dw) CALL histdef(hist_id, 'soiltype', 'Fraction of soil textures', '%', & & iim,jjm, hori_id, nbsoil, 1, nbsoil, soltax_id, 32, once(8), dt,dw) !- !- SECHIBA_HISTLEVEL = 10 !- CALL histdef(hist_id, 'cimean', 'Stomatal CO2 concentation', 'mmole/m2/s', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw) CALL histdef(hist_id, 'vbeta3', 'Beta for Transpiration', 'mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw) CALL histdef(hist_id, 'rveget', 'Surface resistance', 'mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw) CALL histdef(hist_id,'vbeta2','Beta for Interception loss','mm/d', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw) CALL histdef(hist_id, 'qsintmax', 'Maximum Interception capacity', 'Kg/m^2', & & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw) !- CALL histend(hist_id) !===================================================================== !- 3.2 STOMATE's history file !===================================================================== IF ( control_flags%ok_co2 ) THEN !- ! STOMATE IS ACTIVATED !- !Config Key = STOMATE_OUTPUT_FILE !Config Desc = Name of file in which STOMATE's output is going !Config to be written !Config Def = stomate_history.nc !Config Help = This file is going to be created by the model !Config and will contain the output from the model. !Config This file is a truly COADS compliant netCDF file. !Config It will be generated by the hist software from !Config the IOIPSL package. !- stom_histname='stomate_history.nc' CALL getin('STOMATE_OUTPUT_FILE', stom_histname) WRITE(*,*) 'STOMATE_OUTPUT_FILE', stom_histname !- !Config Key = STOMATE_HIST_DT !Config Desc = STOMATE history time step (d) !Config Def = 10. !Config Help = Time step of the STOMATE history file !- hist_days_stom = 10. CALL getin('STOMATE_HIST_DT', hist_days_stom) hist_dt_stom = NINT( hist_days_stom ) * un_jour WRITE(*,*) 'output frequency for STOMATE history file (d): ', & hist_dt_stom/un_jour !- initialize CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & & istp_old, date0, dt, hori_id, hist_id_stom) !- define PFT axis hist_PFTaxis = (/ ( FLOAT(i), i=1,npft ) /) !- declare this axis CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', & & '-', npft, hist_PFTaxis, hist_PFTaxis_id) !- define STOMATE history file CALL stom_define_history (hist_id_stom, npft, iim, jjm, & & dt, hist_dt_stom, hori_id, hist_PFTaxis_id) !- end definition CALL histend(hist_id_stom) !- ENDIF RETURN END SUBROUTINE intsurf_history SUBROUTINE stom_define_history & & (hist_id_stom, npft, iim, jjm, dt, & & hist_dt, hist_hori_id, hist_PFTaxis_id) !--------------------------------------------------------------------- !- Tell ioipsl which variables are to be written !- and on which grid they are defined !--------------------------------------------------------------------- !- !- Input !- !- File id INTEGER,INTENT(in) :: hist_id_stom !- number of PFTs INTEGER,INTENT(in) :: npft !- Domain size INTEGER,INTENT(in) :: iim, jjm !- Time step of STOMATE (seconds) REAL,INTENT(in) :: dt !- Time step of history file (s) REAL,INTENT(in) :: hist_dt !- id horizontal grid INTEGER,INTENT(in) :: hist_hori_id !- id of PFT axis INTEGER,INTENT(in) :: hist_PFTaxis_id !- !- 1 local !- !- maximum history level INTEGER, PARAMETER :: max_hist_level = 10 !- output level (between 0 and 10) !- ( 0:nothing is written, 10:everything is written) INTEGER :: hist_level !- Character strings to define operations for histdef CHARACTER(LEN=40),DIMENSION(max_hist_level) :: ave !- !- Define variables to histdef !- INTEGER :: i,k1,k2,k3,k4,k5 !- Number of variables INTEGER,PARAMETER :: n_v = 68 !- Names of variables CHARACTER(LEN=20),DIMENSION(n_v) :: v_n = (/ & 'SPACE_NAT ','LITTER_STR_AB_NAT ','LITTER_STR_AB_AGRI ',& 'LITTER_MET_AB_NAT ','LITTER_MET_AB_AGRI ','LITTER_STR_BE_NAT ',& 'LITTER_STR_BE_AGRI ','LITTER_MET_BE_NAT ','LITTER_MET_BE_AGRI ',& 'DEADLEAF_COVER ','CARBON_ACTIVE_NAT ','CARBON_ACTIVE_AGRI ',& 'CARBON_SLOW_NAT ','CARBON_SLOW_AGRI ','CARBON_PASSIVE_NAT ',& 'CARBON_PASSIVE_AGRI ','T2M_MONTH ','T2M_WEEK ',& 'HET_RESP_NAT ','HET_RESP_AGRI ','BLACK_CARBON ',& 'FIREFRAC_NAT ','FIREFRAC_AGRI ','FIREINDEX_NAT ',& 'LITTERHUM ','LAI ','VEGET ',& 'VEGET_MAX ','NPP ','GPP ',& 'IND ','LEAF_M ','SAP_M_AB ',& 'SAP_M_BE ','HEART_M_AB ','HEART_M_BE ',& 'ROOT_M ','FRUIT_M ','RESERVE_M ',& 'LEAF_TURN ','SAP_AB_TURN ','ROOT_TURN ',& 'FRUIT_TURN ','LEAF_BM_LITTER ','SAP_AB_BM_LITTER ',& 'SAP_BE_BM_LITTER ','HEART_AB_BM_LITTER ','HEART_BE_BM_LITTER ',& 'ROOT_BM_LITTER ','FRUIT_BM_LITTER ','RESERVE_BM_LITTER ',& 'MAINT_RESP ','GROWTH_RESP ','AGE ',& 'HEIGHT ','MOISTRESS ','VCMAX ',& 'LEAF_AGE ','MORTALITY ','FIREDEATH ',& 'IND_ESTAB ','LIGHT_DEATH ','BM_ALLOC_LEAF ',& 'BM_ALLOC_SAP_AB ','BM_ALLOC_SAP_BE ','BM_ALLOC_ROOT ',& 'BM_ALLOC_FRUIT ','BM_ALLOC_RES ' /) !--------------------------------- !- Descriptions of variables CHARACTER(LEN=50),DIMENSION(n_v) :: v_d = (/ & 'fraction of total space that is natural ',& 'structural litter above nat. ground ',& 'structural litter above agric. ground ',& 'metabolic litter above nat. ground ',& 'metabolic litter above agric. ground ',& 'structural litter below nat. ground ',& 'structural litter below agric. ground ',& 'metabolic litter below nat. ground ',& 'metabolic litter below agric. ground ',& 'fraction of soil covered by dead leaves ',& 'active soil carbon in nat. ground ',& 'active soil carbon in agric. ground ',& 'slow soil carbon in nat. ground ',& 'slow soil carbon in agric. ground ',& 'passive soil carbon in nat. ground ',& 'passive soil carbon in agric. ground ',& 'Monthly 2 m temperature ',& 'Weekly 2 m temperature ',& 'heterotr. resp. from nat. ground ',& 'heterotr. resp. from agric. ground ',& 'black carbon on average total ground ',& 'Fire fraction on natural ground ',& 'Fire fraction on agricultural ground ',& 'Fire index on natural ground ',& 'litter humidity ',& 'Leaf Area Index ',& 'Vegetation fraction ',& 'Maximum vegetation fraction (LAI -> infinity) ',& 'Net primary productivity ',& 'Gross primary productivity ',& 'Density of individuals ',& 'Leaf mass ',& 'Sap mass above ground ',& 'Sap mass below ground ',& 'Heartwood mass above ground ',& 'Heartwood mass below ground ',& 'Root mass ',& 'Fruit mass ',& 'Carbohydrate reserve mass ',& 'Leaf turnover ',& 'Sap turnover above ',& 'Root turnover ',& 'Fruit turnover ',& 'Leaf death ',& 'Sap death above ground ',& 'Sap death below ground ',& 'Heartwood death above ground ',& 'Heartwood death below ground ',& 'Root death ',& 'Fruit death ',& 'Carbohydrate reserve death ',& 'Maintenance respiration ',& 'Growth respiration ',& 'age ',& 'height ',& 'weekly moisture stress ',& 'Maximum rate of carboxylation ',& 'leaf age ',& 'Fraction of trees that dies (gap) ',& 'Fraction of plants killed by fire ',& 'Density of newly established saplings ',& 'Fraction of plants that dies (light competition) ',& 'biomass allocated to leaves ',& 'biomass allocated to sapwood above ground ',& 'biomass allocated to sapwood below ground ',& 'biomass allocated to roots ',& 'biomass allocated to fruits ',& 'biomass allocated to carbohydrate reserve ' /) !- Units of variables CHARACTER(LEN=20),DIMENSION(n_v) :: v_u = (/ & '- ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& 'gC/m^2 (n/a) ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& 'gC/m^2 (n/a) ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& '- ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& 'gC/m^2 (n/a) ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& 'gC/m^2 (n/a) ','K ','K ',& 'gC/m^2 tot/day ','gC/m^2 tot/day ','gC/m^2 tot ',& '1/day ','1/day ','- ',& '- ','- ','- ',& '- ','gC/day/(m^2 (n/a)) ','gC/day/(m^2 (n/a)) ',& '1/(m^2 (n/a)) ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& 'gC/m^2 (n/a) ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& 'gC/m^2 (n/a) ','gC/m^2 (n/a) ','gC/m^2 (n/a) ',& 'gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ',& 'gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ',& 'gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ',& 'gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ','gC/m^2 (n/a)/day ',& 'gC/m^2 tot/day ','gC/m^2 tot/day ','years ',& 'm ','- ','- ',& 'days ','1/day ','1/day ',& '1/day ','1/day ','gC/m**2 nat/agri/dt ',& 'gC/m**2 nat/agri/dt ','gC/m**2 nat/agri/dt ','gC/m**2 nat/agri/dt ',& 'gC/m**2 nat/agri/dt ','gC/m**2 nat/agri/dt ' /) !- Dimensions of variables INTEGER,DIMENSION(n_v) :: v_t = (/ (2,i=1,25), (3,i=26,n_v) /) !- Level of operations for variables INTEGER,DIMENSION(n_v) :: v_l = (/ & 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, & 5, 5, 5, 5, 5, 5, 1, 1, 3, 3, & 10, 1, 1, 10, 5, 1, 1, 1, 1, 1, & 3, 2, 2, 2, 2, 2, 2, 2, 2, 4, & 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & 4, 2, 2, 7, 7, 3, 6, 6, 6, 6, & 6, 6, 5, 5, 5, 5, 5, 5 /) !--------------------------------------------------------------------- !===================================================================== !- 1 history level !===================================================================== !- 1.1 define history level !===================================================================== !Config Key = STOMATE_HISTLEVEL !Config Desc = STOMATE history output level (0..10) !Config Def = 10 !Config Help = 0: nothing is written; 10: everything is written !- hist_level = 10 CALL getin('STOMATE_HISTLEVEL', hist_level) !- WRITE(*,*) 'STOMATE history level: ',hist_level IF ( (hist_level > max_hist_level).OR.(hist_level < 0) ) THEN STOP 'This history level is not allowed' ENDIF !===================================================================== !- 1.2 define operations according to output level !===================================================================== ave(1:10) = (/ ('ave(scatter(X))',i=1,hist_level), & & ('never ',i=hist_level+1,10) /) !===================================================================== !- 2 surface fields (2d) !- 3 PFT: 3rd dimension !===================================================================== k2 = 1; k5 = 32; DO i=1,n_v IF (v_t(i) == 2) THEN k1 = 1; k3 = 1; k4 = -99; ELSE IF (v_t(i) == 3) THEN k1 = npft; k3 = npft; k4 = hist_PFTaxis_id; ELSE CYCLE ENDIF CALL histdef (hist_id_stom, & & TRIM(v_n(i)),TRIM(v_d(i)),TRIM(v_u(i)), iim,jjm, & & hist_hori_id, k1,k2,k3,k4,k5, ave(v_l(i)), dt, hist_dt) ENDDO !--------------------------------- END SUBROUTINE stom_define_history SUBROUTINE grid_stuff (npts, iim, jjm, grid_lon, grid_lat, kindex, & neighbours, resolution) ! ! 0 interface ! ! ! 0.1 input ! ! Domain size INTEGER, INTENT(in) :: npts ! Size of cartesian grid INTEGER, INTENT(in) :: iim, jjm ! Longitudes on cartesian grid REAL, DIMENSION(iim,jjm), INTENT(in) :: grid_lon ! Latitudes on cartesian grid REAL, DIMENSION(iim,jjm), INTENT(in) :: grid_lat ! Indices of STOMATE points on cartesian grid INTEGER, DIMENSION(npts), INTENT(in) :: kindex ! ! 0.2 output ! ! indices of the 4 neighbours of each grid point (1=N, 2=E, 3=S, 4=W) ! a zero or negative index means that this neighbour is not a land point INTEGER, DIMENSION(npts,4), INTENT(out) :: neighbours ! resolution at each grid point in m (1=E-W, 2=N-S) REAL, DIMENSION(npts,2), INTENT(out) :: resolution ! ! 0.3 local ! ! default resolution (m) REAL :: default_resolution ! which STOMATE point corresponds to the given point on the cartesian grid INTEGER, DIMENSION(iim,jjm) :: correspondance ! cosine of the latitude REAL :: coslat ! number of points where default resolution is used INTEGER :: ndefault_lon, ndefault_lat ! Indices INTEGER :: i,ip,jp REAL, PARAMETER :: R_Earth = 6378000. INTEGER, SAVE :: bavard=1 REAL, PARAMETER :: min_stomate = 1.E-8 REAL, SAVE :: pi ! ========================================================================= pi = 4. * ATAN(1.) IF ( bavard .GE. 4 ) WRITE(*,*) 'Entering grid_stuff' ! default resolution default_resolution = 250000. IF ( bavard .GT. 1 ) WRITE(*,*) 'grid stuff: default resolution (m): ',default_resolution ndefault_lon = 0 ndefault_lat = 0 ! initialize output neighbours(:,:) = -1 resolution(:,:) = 0. correspondance(:,:) = -1 DO i = 1, npts ! ! 1 find numbers of the latitude and longitude of each point ! ! index of latitude jp = INT( (kindex(i)-1) /iim ) + 1 ! index of longitude ip = kindex(i) - ( jp-1 ) * iim correspondance(ip,jp) = kindex(i) ENDDO DO i = 1, npts ! index of latitude jp = INT( (kindex(i)-1) /iim ) + 1 ! index of longitude ip = kindex(i) - ( jp-1 ) * iim ! ! 2 resolution ! ! ! 2.1 longitude ! ! prevent infinite resolution at the pole coslat = AMAX1( COS( grid_lat(ip,jp) * pi/180. ), 0.001 ) IF ( iim .GT. 1 ) THEN IF ( ip .EQ. 1 ) THEN resolution(i,1) = ABS( grid_lon(ip+1,jp) - grid_lon(ip,jp) ) * & pi/180. * R_Earth * coslat ELSEIF ( ip .EQ. iim ) THEN resolution(i,1) = ABS( grid_lon(ip,jp) - grid_lon(ip-1,jp) ) * & pi/180. * R_Earth * coslat ELSE resolution(i,1) = ABS( grid_lon(ip+1,jp) - grid_lon(ip-1,jp) ) / 2. * & pi/180. * R_Earth * coslat ENDIF ELSE resolution(i,1) = default_resolution ndefault_lon = ndefault_lon + 1 ENDIF ! ! 2.2 latitude ! IF ( jjm .GT. 1 ) THEN IF ( jp .EQ. 1 ) THEN resolution(i,2) = ABS( grid_lat(ip,jp) - grid_lat(ip,jp+1) ) * pi/180. * R_Earth ELSEIF ( jp .EQ. jjm ) THEN resolution(i,2) = ABS( grid_lat(ip,jp-1) - grid_lat(ip,jp) ) * pi/180. * R_Earth ELSE resolution(i,2) = ABS( grid_lat(ip,jp-1) - grid_lat(ip,jp+1) ) / 2. * pi/180. * R_Earth ENDIF ELSE resolution(i,2) = default_resolution ndefault_lat = ndefault_lat + 1 ENDIF ! ! 3 find neighbours ! ! North IF ( jp .GT. 1 ) THEN neighbours(i,1) = correspondance(ip,jp-1) ELSE ! special treatment for the pole if we are really in a 2d grid IF ( ( iim .GT. 1 ) .AND. ( jjm .GT. 1 ) ) THEN IF ( ABS( ( grid_lat(ip,jp) ) - 90. ) .LT. min_stomate ) THEN ! the grid point sits exactly on the pole. The neighbour is situated ! at a lower latitude. neighbours(i,1) = correspondance( MOD(ip+iim/2-1,iim)+1, jp+1 ) ELSE ! look across the North Pole neighbours(i,1) = correspondance( MOD(ip+iim/2-1,iim)+1, jp ) ENDIF ENDIF ENDIF ! South IF ( jp .LT. jjm ) THEN neighbours(i,3) = correspondance(ip,jp+1) ELSE ! special treatment for the pole if we are really in a 2d grid IF ( ( iim .GT. 1 ) .AND. ( jjm .GT. 1 ) ) THEN IF ( ( ABS( grid_lat(ip,jp) ) - 90. ) .LT. min_stomate ) THEN ! the grid point sits exactly on the pole. The neighbour is situated ! at a lower latitude. neighbours(i,3) = correspondance( MOD(ip+iim/2-1,iim)+1, jp-1 ) ELSE ! look across the South Pole neighbours(i,3) = correspondance( MOD(ip+iim/2-1,iim)+1, jp ) ENDIF ENDIF ENDIF ! East IF ( ip .LT. iim ) THEN neighbours(i,2) = correspondance(ip+1,jp) ELSEIF ( iim .GT. 1 ) THEN neighbours(i,2) = correspondance(1,jp) ENDIF ! West IF ( ip .GT. 1 ) THEN neighbours(i,4) = correspondance(ip-1,jp) ELSEIF ( iim .GT. 1 ) THEN neighbours(i,4) = correspondance(iim,jp) ENDIF ENDDO IF ( bavard .GT. 1 ) THEN WRITE(*,*) ' > Total number of points: ',npts WRITE(*,*) ' > Using default zonal resolution at',ndefault_lon,' points.' WRITE(*,*) ' > Using default meridional resolution at',ndefault_lat,' points.' ENDIF IF ( bavard .GT. 1 ) WRITE(*,*) 'Leaving grid_stuff' END SUBROUTINE grid_stuff END MODULE intersurf