source: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 @ 3203

Last change on this file since 3203 was 3181, checked in by jbclement, 10 months ago

PEM:

  • Addition of a script "inipem_orbit.sh" in the deftank to modify the orbital parameters of a file "startfi.nc" according to the date set in the file "run_PEM.def" and data found in "obl_ecc_lsp.asc";
  • Flow of glaciers is now computed only when there are slopes;
  • Reversion to the name "diagpem.nc" for the PEM outputs (as decided) which was modified in r3171;
  • Some small cleanings.

JBC

File size: 15.6 KB
RevLine 
[3096]1MODULE read_data_PCM_mod
[2794]2
[3076]3use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
4                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
5                  nf90_inquire_dimension, nf90_close
[2794]6
[3076]7implicit none
8
[3122]9character(13), parameter :: modname = 'read_data_PCM'
10character(256)           :: msg      ! for reading
11integer                  :: fID, vID ! for reading
[3076]12
[2794]13!=======================================================================
[3076]14contains
15!=======================================================================
16
[3149]17SUBROUTINE read_data_PCM(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM_phys,ps_timeseries,          &
18                         min_co2_ice,min_h2o_ice,tsurf_avg,tsoil_avg,tsurf_PCM,tsoil_PCM,q_co2,q_h2o,co2_ice_slope, &
19                         watersurf_density_avg,watersoil_density)
[3076]20use comsoil_h,             only: nsoilmx
21use comsoil_h_PEM,         only: soil_pem
22use constants_marspem_mod, only: m_co2, m_noco2
23
24implicit none
25
26!=======================================================================
[2794]27!
[3096]28! Purpose: Read initial confitions file from the PCM
[2794]29!
[2855]30! Authors: RV & LL
[2794]31!=======================================================================
32
[3076]33include "dimensions.h"
[2794]34
[3076]35!=======================================================================
[2794]36! Arguments:
[3149]37character(*), intent(in) :: filename                            ! File name
38integer,      intent(in) :: timelen                             ! number of times stored in the file
39integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
[2855]40! Ouputs
[3181]41real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice per slope of the year [kg/m^2]
[3076]42real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
[3181]43real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_PCM_phys ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
[3076]44real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
45real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
46real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
[3149]47real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per slope of the year [kg/m^2]
[2794]48!SOIL
[3149]49real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg             ! Average surface temperature of the concatenated file [K]
50real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg             ! Average soil temperature of the concatenated file [K]
51real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_PCM             ! Surface temperature of the concatenated file, time series [K]
52real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_PCM             ! Soil temperature of the concatenated file, time series [K]
53real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg ! Water density at the surface [kg/m^3]
[3076]54real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
[3122]55!=======================================================================
56! Local Variables
[3076]57character(12)                   :: start_file_type = "earth" ! default start file type
58real, dimension(:), allocatable :: time      ! times stored in start
59integer                         :: indextime ! index of selected time
60integer                         :: edges(4), corner(4)
61integer                         :: i, j, l, t          ! loop variables
[3122]62real                            :: A, B, mmean         ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
[3076]63integer                         :: islope              ! loop for variables
64character(2)                    :: num                 ! for reading sloped variables
[3123]65real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! h2o ice per slope of the concatenated file [kg/m^2]
[3122]66real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
[3130]67real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
[3149]68real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_PCM               ! CO2 volume mixing ratio in the first layer [mol/m^3]
69real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_PCM                    ! Surface Pressure [Pa]
[3076]70real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
71real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
[3149]72real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_avg_dyn             ! Average surface temperature of the concatenated file [K]
73real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Average soil temperature of the concatenated file [K]
74real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_PCM_dyn             ! Surface temperature of the concatenated file, time series [K]
75real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_PCM_dyn             ! Soil temperature of the concatenated file, time series [K]
[3123]76real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
77real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
78real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn         ! co2 ice amount per  slope of the year [kg/m^2]
79real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn     ! Water density at the surface, time series [kg/m^3]
[3149]80real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_avg ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
[3123]81real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
82real, dimension(ngrid,nslope,timelen)                               :: watersurf_density         ! Water density at the surface, physical grid, time series [kg/m^3]
[2794]83!-----------------------------------------------------------------------
[3122]84! Open the NetCDF file
85write(*,*) "Opening "//filename//"..."
86call error_msg(NF90_OPEN(filename,NF90_NOWRITE,fID),"open",filename)
[2794]87
[3122]88! Dowload the data from the file
[3149]89call get_var3("ps",ps_PCM)
[3122]90write(*,*) "Data for surface pressure downloaded!"
[2794]91
[3122]92call get_var3("co2_layer1",q_co2_dyn)
93write(*,*) "Data for vmr co2 downloaded!"
[2794]94
[3122]95call get_var3("h2o_layer1",q_h2o_dyn)
96write(*,*) "Data for vmr h2o downloaded!"
[2794]97
[3122]98if (nslope == 1) then ! There is no slope
99    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
100    write(*,*) "Data for co2_ice downloaded!"
[2835]101
[3122]102    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
103    write(*,*) "Data for h2o_ice_s downloaded!"
[2835]104
[3149]105    call get_var3("tsurf",tsurf_PCM_dyn(:,:,1,:))
[3143]106    write(*,*) "Data for tsurf downloaded!"
107
[3122]108#ifndef CPP_STD
109    call get_var3("watercap",watercap(:,:,1,:))
110    write(*,*) "Data for watercap downloaded!"
[3130]111
112    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
113    write(*,*) "Data for perennial_co2ice downloaded!"
[2835]114
[3122]115    if (soil_pem) then
[3149]116        call get_var4("soiltemp",tsoil_PCM_dyn(:,:,:,1,:))
[3122]117        write(*,*) "Data for soiltemp downloaded!"
[2794]118
[3152]119        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
[3122]120        write(*,*) "Data for waterdensity_soil downloaded!"
[2794]121
[3152]122        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
[3122]123        write(*,*) "Data for waterdensity_surface downloaded!"
124    endif !soil_pem
125#endif
126else ! We use slopes
[3076]127    do islope = 1,nslope
128        write(num,'(i2.2)') islope
129        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
130    enddo
[3122]131    write(*,*) "Data for co2_ice downloaded!"
132
[3076]133    do islope = 1,nslope
134        write(num,'(i2.2)') islope
135        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
136    enddo
[3122]137    write(*,*) "Data for h2o_ice_s downloaded!"
[2980]138
[3143]139    do islope = 1,nslope
140        write(num,'(i2.2)') islope
[3149]141        call get_var3("tsurf_slope"//num,tsurf_PCM_dyn(:,:,islope,:))
[3143]142    enddo
143    write(*,*) "Data for tsurf downloaded!"
144
[2985]145#ifndef CPP_STD
[3076]146    do islope = 1,nslope
147        write(num,'(i2.2)') islope
[3122]148        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
[3076]149    enddo
[3122]150    write(*,*) "Data for watercap downloaded!"
[3130]151
152    do islope = 1,nslope
153        write(num,'(i2.2)') islope
154        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
155    enddo
156    write(*,*) "Data for perennial_co2ice downloaded!"
[2985]157
[3076]158    if (soil_pem) then
159        do islope = 1,nslope
160            write(num,'(i2.2)') islope
[3149]161            call get_var4("soiltemp_slope"//num,tsoil_PCM_dyn(:,:,:,islope,:))
[3076]162        enddo
[3122]163        write(*,*) "Data for soiltemp downloaded!"
[2835]164
[3076]165        do islope = 1,nslope
166            write(num,'(i2.2)') islope
[3122]167            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
[3076]168        enddo
[3122]169        write(*,*) "Data for waterdensity_soil downloaded!"
[2794]170
[3076]171        do islope = 1,nslope
172            write(num,'(i2.2)') islope
[3122]173            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
[3076]174        enddo
[3122]175        write(*,*) "Data for waterdensity_surface downloaded!"
[3076]176    endif !soil_pem
[2985]177#endif
[3122]178endif
[2897]179
[2794]180! Compute the minimum over the year for each point
[3122]181write(*,*) "Computing the min of h2o_ice..."
182where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
[3149]183min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
[3122]184write(*,*) "Computing the min of co2_ice..."
185where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
[3149]186min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
[2794]187
[3122]188! Compute averages over the year for each point
189write(*,*) "Computing the average of tsurf..."
[3149]190tsurf_avg_dyn = sum(tsurf_PCM_dyn,4)/timelen
[2794]191
[3076]192if (soil_pem) then
[3122]193    write(*,*) "Computing average of tsoil..."
[3149]194    tsoil_avg_dyn = sum(tsoil_PCM_dyn,5)/timelen
[3122]195    write(*,*) "Computing average of waterdensity_surface..."
[3149]196    watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
[3076]197endif
[2794]198
[3122]199! By definition, we get rid of the negative values
200A = (1./m_co2 - 1./m_noco2)
201B = 1./m_noco2
[3076]202do i = 1,iim + 1
203    do j = 1,jjm_input + 1
204        do t = 1, timelen
205            if (q_co2_dyn(i,j,t) < 0) then
206                q_co2_dyn(i,j,t) = 1.e-10
207            else if (q_co2_dyn(i,j,t) > 1) then
208                q_co2_dyn(i,j,t) = 1.
209            endif
210            if (q_h2o_dyn(i,j,t) < 0) then
[3130]211                q_h2o_dyn(i,j,t) = 1.e-10
[3076]212            else if (q_h2o_dyn(i,j,t) > 1) then
213                q_h2o_dyn(i,j,t) = 1.
214            endif
215            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
[3149]216            vmr_co2_PCM(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
[3076]217        enddo
218    enddo
219enddo
[2794]220
[2980]221#ifndef CPP_1D
[3149]222    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_PCM,vmr_co2_PCM_phys)
223    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_PCM,ps_timeseries)
[3076]224    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
225    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
[2980]226
[3076]227    do islope = 1,nslope
228        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
229        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
230        if (soil_pem) then
231            do l = 1,nsoilmx
[3149]232                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
[3122]233                do t = 1,timelen
[3149]234                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_PCM_dyn(:,:,l,islope,t),tsoil_PCM(:,l,islope,t))
[3076]235                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
236                enddo
237            enddo
[3149]238            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
[3122]239        endif ! soil_pem
[3076]240        do t = 1,timelen
[3149]241            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_PCM_dyn(:,:,islope,t),tsurf_PCM(:,islope,t))
[3076]242            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
243        enddo
244    enddo
[2849]245
[3149]246    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
[3076]247#else
[3149]248    vmr_co2_PCM_phys(1,:) = vmr_co2_PCM(1,1,:)
249    ps_timeseries(1,:) = ps_PCM(1,1,:)
[3076]250    q_co2(1,:) = q_co2_dyn(1,1,:)
251    q_h2o(1,:) = q_h2o_dyn(1,1,:)
252    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
253    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
254    if (soil_pem) then
[3149]255        tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
256        tsoil_PCM(1,:,:,:) = tsoil_PCM_dyn(1,1,:,:,:)
[3076]257        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
[3149]258        watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
[3122]259    endif ! soil_pem
[3149]260    tsurf_PCM(1,:,:) = tsurf_PCM_dyn(1,1,:,:)
[3076]261    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
[3149]262    tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
[3076]263#endif
[2897]264
[3096]265END SUBROUTINE read_data_PCM
[2897]266
[3122]267!=======================================================================
268
[3076]269SUBROUTINE check_dim(n1,n2,str1,str2)
[2980]270
[3076]271implicit none
[2980]272
[3076]273integer,            intent(in) :: n1, n2
274character(len = *), intent(in) :: str1, str2
[2980]275
[3076]276character(256) :: s1, s2
[2794]277
[3076]278if (n1 /= n2) then
279    s1 = 'value of '//trim(str1)//' ='
280    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
281    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
282    call abort_gcm(trim(modname),trim(msg),1)
283endif
284
[2794]285END SUBROUTINE check_dim
286
[3076]287!=======================================================================
[2794]288
[3130]289SUBROUTINE get_var1(var,v)
[3076]290
291implicit none
292
[3130]293character(len = *), intent(in)  :: var
[3076]294real, dimension(:), intent(out) :: v
295
[3130]296call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
297call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
[3076]298
[2794]299END SUBROUTINE get_var1
300
[3076]301!=======================================================================
[2794]302
[3130]303SUBROUTINE get_var3(var,v) ! on U grid
[2794]304
[3076]305implicit none
306
[3130]307character(len = *),     intent(in)  :: var
[3076]308real, dimension(:,:,:), intent(out) :: v
309
[3130]310call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
311call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
[3076]312
[2794]313END SUBROUTINE get_var3
314
[3076]315!=======================================================================
316
[3130]317SUBROUTINE get_var4(var,v)
[3076]318
319implicit none
320
[3130]321character(len = *),       intent(in)  :: var
[3076]322real, dimension(:,:,:,:), intent(out) :: v
323
[3130]324call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
325call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
[3076]326
[2794]327END SUBROUTINE get_var4
328
[3076]329!=======================================================================
[2794]330
[3076]331SUBROUTINE error_msg(ierr,typ,nam)
332
333implicit none
334
335integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
336character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
337character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
338
339if (ierr == nf90_noerr) return
340select case(typ)
341    case('inq');   msg="Field <"//trim(nam)//"> is missing"
342    case('get');   msg="Reading failed for <"//trim(nam)//">"
343    case('open');  msg="File opening failed for <"//trim(nam)//">"
344    case('close'); msg="File closing failed for <"//trim(nam)//">"
345    case default
346        write(*,*) 'There is no message for this error.'
347        error stop
348end select
349call abort_gcm(trim(modname),trim(msg),ierr)
350
351END SUBROUTINE error_msg
352
[3096]353END MODULE read_data_PCM_mod
Note: See TracBrowser for help on using the repository browser.