! ! $Id $ ! SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,min_h2o_ice_s,min_co2_ice_s,vmr_co2_gcm,ps_GCM, & min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope, & watersurf_density,watersoil_density) use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, & nf90_get_var, nf90_inq_varid, nf90_inq_dimid, & nf90_inquire_dimension,nf90_close use comsoil_h, only: nsoilmx USE comsoil_h_PEM, ONLY: soil_pem IMPLICIT NONE !======================================================================= ! ! Purpose: Read initial confitions file from the GCM ! ! Authors: RV & LL !======================================================================= include "dimensions.h" !=============================================================================== ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: fichnom !--- FILE NAME INTEGER, INTENT(IN) :: timelen ! number of times stored in the file INTEGER :: iim_input,jjm_input,nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes ! Ouputs REAL, INTENT(OUT) :: min_h2o_ice_s(iim_input+1,jjm_input+1) ! Minimum of h2o ice, mesh averaged of the year [kg/m^2] REAL, INTENT(OUT) :: min_co2_ice_s(iim_input+1,jjm_input+1) ! Minimum of co2 ice, mesh averaged of the year [kg/m^2] REAL, INTENT(OUT) :: min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2 ice per slope of the year [kg/m^2] REAL, INTENT(OUT) :: min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2] REAL, INTENT(OUT) :: vmr_co2_gcm(iim_input+1,jjm_input+1,timelen) ! CO2 volume mixing ratio in the first layer [mol/m^3] REAL, INTENT(OUT) :: q_h2o_GCM(iim_input+1,jjm_input+1,timelen) ! H2O mass mixing ratio in the first layer [kg/m^3] REAL, INTENT(OUT) :: q_co2_GCM(iim_input+1,jjm_input+1,timelen) ! CO2 mass mixing ratio in the first layer [kg/m^3] REAL, INTENT(OUT) :: ps_GCM(iim_input+1,jjm_input+1,timelen) ! Surface Pressure [Pa] REAL, INTENT(OUT) :: co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! co2 ice amount per slope of the year [kg/m^2] !SOIL REAL, INTENT(OUT) :: tsurf_ave(iim_input+1,jjm_input+1,nslope) ! Average surface temperature of the concatenated file [K] REAL, INTENT(OUT) :: tsoil_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average soil temperature of the concatenated file [K] REAL ,INTENT(OUT) :: tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file, time series [K] REAL , INTENT(OUT) :: tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file, time series [K] REAL , INTENT(OUT) :: watersurf_density(iim_input+1,jjm_input+1,nslope,timelen) ! Water density at the surface, time series [kg/m^3] REAL , INTENT(OUT) :: watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3] REAL :: TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia of the concatenated file, times series [SI] REAL, INTENT(OUT) :: TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average Thermal Inertia of the concatenated file [SI] !=============================================================================== ! Local Variables CHARACTER(LEN=256) :: msg, var, modname ! for reading INTEGER :: iq, fID, vID, idecal ! for reading INTEGER :: ierr ! for reading CHARACTER(len=12) :: start_file_type="earth" ! default start file type REAL,ALLOCATABLE :: time(:) ! times stored in start INTEGER :: indextime ! index of selected time INTEGER :: edges(4),corner(4) INTEGER :: i,j,t ! loop variables real,save :: m_co2, m_noco2, A , B, mmean ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer INTEGER :: islope ! loop for variables CHARACTER*2 :: num ! for reading sloped variables REAL, ALLOCATABLE :: h2o_ice_s(:,:,:) ! h2o ice, mesh averaged, of the concatenated file [kg/m^2] REAL, ALLOCATABLE :: co2_ice_s(:,:,:) ! co2 ice, mesh averaged, of the concatenated file [kg/m^2] REAL, ALLOCATABLE :: h2o_ice_s_slope(:,:,:,:) ! h2o ice per slope of the concatenated file [kg/m^2] REAL, ALLOCATABLE :: watercap_slope(:,:,:,:) !----------------------------------------------------------------------- modname="read_data_gcm" m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) A =(1/m_co2 - 1/m_noco2) B=1/m_noco2 allocate(co2_ice_s(iim+1,jjm+1,timelen)) allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen)) allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen)) allocate(h2o_ice_s(iim+1,jjm+1,timelen)) allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen)) print *, "Opening ", fichnom, "..." ! Open initial state NetCDF file var=fichnom CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var) print *, "Downloading data for h2oice ..." ! Get h2o_ice_s of the concatenated file CALL get_var3("h2o_ice_s" ,h2o_ice_s) print *, "Downloading data for h2oice done" print *, "Downloading data for co2ice ..." CALL get_var3("co2ice" ,co2_ice_s) print *, "Downloading data for co2ice done" print *, "Downloading data for vmr co2..." CALL get_var3("co2_cropped" ,q_co2_GCM) print *, "Downloading data for vmr co2 done" print *, "Downloading data for vmr h20..." CALL get_var3("h2o_cropped" ,q_h2o_GCM) print *, "Downloading data for vmr h2o done" print *, "Downloading data for surface pressure ..." CALL get_var3("ps" ,ps_GCM) print *, "Downloading data for surface pressure done" print *, "nslope=", nslope print *, "Downloading data for co2ice_slope ..." if(nslope.gt.1) then DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var3("co2ice_slope"//num,co2_ice_slope(:,:,islope,:)) ENDDO print *, "Downloading data for co2ice_slope done" print *, "Downloading data for h2o_ice_s_slope ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_slope(:,:,islope,:)) ENDDO print *, "Downloading data for h2o_ice_s_slope done" print *, "Downloading data for watercap_slope ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:)) ENDDO print *, "Downloading data for watercap_slope done" print *, "Downloading data for tsurf_slope ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var3("tsurf_slope"//num,tsurf_gcm(:,:,islope,:)) ENDDO print *, "Downloading data for tsurf_slope done" if(soil_pem) then print *, "Downloading data for tsoil_slope ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var4("tsoil_slope"//num,tsoil_gcm(:,:,:,islope,:)) ENDDO print *, "Downloading data for tsoil_slope done" print *, "Downloading data for inertiesoil_slope ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var4("inertiesoil_slope"//num,TI_gcm(:,:,:,islope,:)) ENDDO print *, "Downloading data for inertiesoil_slope done" print *, "Downloading data for watersoil_density ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var4("Waterdensity_soil_slope"//num,watersoil_density(:,:,:,islope,:)) ENDDO print *, "Downloading data for watersoil_density done" print *, "Downloading data for watersurf_density ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_var3("Waterdensity_surface"//num,watersurf_density(:,:,islope,:)) ENDDO print *, "Downloading data for watersurf_density done" endif !soil_pem else !nslope=1 no slope, we copy all the values co2_ice_slope(:,:,1,:)=co2_ice_s(:,:,:) h2o_ice_s_slope(:,:,1,:)=h2o_ice_s(:,:,:) call get_var3("tsurf",tsurf_gcm(:,:,1,:)) if(soil_pem) then call get_var4("tsoil",tsoil_gcm(:,:,:,1,:)) call get_var4("inertiesoil",TI_gcm(:,:,:,1,:)) endif !soil_pem endif !nslope=1 ! Compute the minimum over the year for each point print *, "Computing the min of h2o_ice" min_h2o_ice_s(:,:)=minval(h2o_ice_s,3) print *, "Computing the min of co2_ice" min_co2_ice_s(:,:)=minval(co2_ice_s,3) print *, "Computing the min of h2o_ice_slope" min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope+watercap_slope,4) print *, "Computing the min of co2_ice_slope" min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4) !Compute averages print *, "Computing average of tsurf" tsurf_ave(:,:,:)=SUM(tsurf_gcm(:,:,:,:),4)/timelen if(soil_pem) then print *, "Computing average of tsoil" tsoil_ave(:,:,:,:)=SUM(tsoil_gcm(:,:,:,:,:),5)/timelen print *, "Computing average of TI" TI_ave(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen endif ! By definition, a density is positive, we get rid of the negative values DO i=1,iim+1 DO j = 1, jjm+1 if (min_co2_ice_s(i,j).LT.0) then min_h2o_ice_s(i,j) = 0. min_co2_ice_s(i,j) = 0. endif DO islope=1,nslope if (min_co2_ice_slope(i,j,islope).LT.0) then min_co2_ice_slope(i,j,islope) = 0. endif ! if (min_h2o_ice_slope(i,j,islope).LT.0) then ! min_h2o_ice_slope(i,j,islope) = 0. ! endif ENDDO ENDDO ENDDO DO i=1,iim+1 DO j = 1, jjm+1 DO t = 1, timelen if (q_co2_GCM(i,j,t).LT.0) then q_co2_GCM(i,j,t)=1E-10 elseif (q_co2_GCM(i,j,t).GT.1) then q_co2_GCM(i,j,t)=1. endif if (q_h2o_GCM(i,j,t).LT.0) then q_h2o_GCM(i,j,t)=1E-30 elseif (q_h2o_GCM(i,j,t).GT.1) then q_h2o_GCM(i,j,t)=1. endif mmean=1/(A*q_co2_GCM(i,j,t) +B) vmr_co2_gcm(i,j,t) = q_co2_GCM(i,j,t)*mmean/m_co2 ENDDO ENDDO ENDDO deallocate(co2_ice_s) deallocate(h2o_ice_s_slope) deallocate(h2o_ice_s) CONTAINS SUBROUTINE check_dim(n1,n2,str1,str2) INTEGER, INTENT(IN) :: n1, n2 CHARACTER(LEN=*), INTENT(IN) :: str1, str2 CHARACTER(LEN=256) :: s1, s2 IF(n1/=n2) THEN s1='value of '//TRIM(str1)//' =' s2=' read in starting file differs from parametrized '//TRIM(str2)//' =' WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2 CALL ABORT_gcm(TRIM(modname),TRIM(msg),1) END IF END SUBROUTINE check_dim SUBROUTINE get_var1(var,v) CHARACTER(LEN=*), INTENT(IN) :: var REAL, INTENT(OUT) :: v(:) CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) CALL err(NF90_GET_VAR(fID,vID,v),"get",var) END SUBROUTINE get_var1 SUBROUTINE get_var3(var,v) ! on U grid CHARACTER(LEN=*), INTENT(IN) :: var REAL, INTENT(OUT) :: v(:,:,:) CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) CALL err(NF90_GET_VAR(fID,vID,v),"get",var) END SUBROUTINE get_var3 SUBROUTINE get_var4(var,v) CHARACTER(LEN=*), INTENT(IN) :: var REAL, INTENT(OUT) :: v(:,:,:,:) CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) CALL err(NF90_GET_VAR(fID,vID,v),"get",var) END SUBROUTINE get_var4 SUBROUTINE err(ierr,typ,nam) INTEGER, INTENT(IN) :: ierr !--- NetCDF ERROR CODE CHARACTER(LEN=*), INTENT(IN) :: typ !--- TYPE OF OPERATION CHARACTER(LEN=*), INTENT(IN) :: nam !--- FIELD/FILE NAME IF(ierr==NF90_NoERR) RETURN SELECT CASE(typ) CASE('inq'); msg="Field <"//TRIM(nam)//"> is missing" CASE('get'); msg="Reading failed for <"//TRIM(nam)//">" CASE('open'); msg="File opening failed for <"//TRIM(nam)//">" CASE('close'); msg="File closing failed for <"//TRIM(nam)//">" END SELECT CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr) END SUBROUTINE err END SUBROUTINE read_data_gcm