Ignore:
Timestamp:
Jan 5, 2024, 3:37:30 PM (12 months ago)
Author:
llange
Message:

PEM
Add the possibily to output the soil fields during a PEM run in a dedicated diagsoilpem

ll

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3170 r3171  
    189189== 04/01/2024 == LL
    190190Fixing a small bug: the subroutine compute_icetable was always called, even if tthe option 'icetable_equilibrium' was set to false in the run_PEM.def. It is now fixed by adding a flag before the call.
     191
     192== 05/01/2024 == LL
     193Add the possibily to output the soil fields during a PEM run in a dedicated diagsoilpem
  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3170 r3171  
    6868use recomp_tend_co2_slope_mod,  only: recomp_tend_co2_slope
    6969use soil_pem_compute_mod,       only: soil_pem_compute
    70 use writediagpem_mod,           only: writediagpem
     70use writediagpem_mod,           only: writediagfipem, writediagsoilpem
    7171
    7272#ifndef CPP_STD
     
    874874!    II_e Outputs
    875875!------------------------
    876     call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_avg_press_new/))
     876    call writediagfipem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_avg_press_new/))
    877877    do islope = 1,nslope
    878878        write(str2(1:2),'(i2.2)') islope
    879         call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
    880         call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
    881         call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
    882         call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
    883         call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
    884         call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
     879        call writediagfipem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
     880        call writediagfipem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
     881        call writediagfipem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
     882        call writediagfipem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
     883        call writediagfipem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
     884        call writediagfipem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
     885        if(icetable_equilibrium) then
     886            call writediagfipem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
     887            call writediagfipem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
     888        endif
     889        if(soil_pem) then
     890            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
     891            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
     892            if (adsorption_pem) then
     893                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
     894                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
     895            endif                       
     896        endif
    885897    enddo
    886898
  • trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90

    r3122 r3171  
    77!=======================================================================
    88
    9 SUBROUTINE writediagpem(ngrid,nom,titre,unite,dim,px)
     9SUBROUTINE writediagfipem(ngrid,nom,titre,unite,dim,px)
    1010
    1111!  Ecriture de variables diagnostiques au choix dans la physique
     
    8686integer                 :: idim, varid
    8787integer                 :: nid
    88 character(*), parameter :: fichnom = "diagpem.nc"
     88character(*), parameter :: fichnom = "diagfipem.nc"
    8989integer, dimension(4)   :: id
    9090integer, dimension(4)   :: edges, corner
     
    578578if (is_master) ierr= NF_CLOSE(nid)
    579579
    580 END SUBROUTINE writediagpem
     580END SUBROUTINE writediagfipem
     581
     582subroutine writediagsoilpem(ngrid,name,title,units,dimpx,px)
     583
     584! Write variable 'name' to NetCDF file 'diagsoil.nc'.
     585! The variable may be 3D (lon,lat,depth) subterranean field,
     586! a 2D (lon,lat) surface field, or a simple scalar (0D variable).
     587!
     588! Calls to 'writediagsoilpem' can originate from anywhere in the program;
     589! An initialisation of variable 'name' is done if it is the first time
     590! that this routine is called with given 'name'; otherwise data is appended
     591! (yielding the sought time series of the variable)
     592
     593! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
     594
     595use comsoil_h, only: nsoilmx, inertiedat
     596use geometry_mod, only: cell_area
     597use time_phylmdz_mod, only: ecritphy, day_step, iphysiq
     598use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
     599use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat
     600use mod_grid_phy_lmdz, only : grid_type, unstructured
     601
     602implicit none
     603
     604include"netcdf.inc"
     605
     606! Arguments:
     607integer,intent(in) :: ngrid ! number of (horizontal) points of physics grid
     608! i.e. ngrid = 2+(jjm-1)*iim - 1/jjm
     609character(len=*),intent(in) :: name ! 'name' of the variable
     610character(len=*),intent(in) :: title ! 'long_name' attribute of the variable
     611character(len=*),intent(in) :: units ! 'units' attribute of the variable
     612integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
     613real,dimension(ngrid,nsoilmx),intent(in) :: px ! variable
     614
     615! Local variables:
     616real*4,dimension(nbp_lon+1,nbp_lat,nsoilmx) :: data3 ! to store 3D data
     617real*4,dimension(nbp_lon+1,nbp_lat) :: data2 ! to store 2D data
     618real*4 :: data0 ! to store 0D data
     619real*4 :: data3_1d(1,nsoilmx) ! to store a profile in 1D model
     620real*4 :: data2_1d ! to store surface value with 1D model
     621integer :: i,j,l ! for loops
     622integer :: ig0
     623
     624real*4,save :: date ! time counter (in elapsed days)
     625
     626real :: inertia((nbp_lon+1),nbp_lat,nsoilmx)
     627real :: area((nbp_lon+1),nbp_lat)
     628
     629real :: inertiafi_glo(klon_glo,nsoilmx)
     630real :: areafi_glo(klon_glo)
     631
     632integer,save :: isample ! sample rate at which data is to be written to output
     633integer,save :: ntime=0 ! counter to internally store time steps
     634character(len=20),save :: firstname="1234567890"
     635integer,save :: zitau=0
     636!$OMP THREADPRIVATE(date,isample,ntime,firstname,zitau)
     637
     638character(len=30) :: filename="diagsoilpem.nc"
     639
     640! NetCDF stuff:
     641integer :: nid ! NetCDF output file ID
     642integer :: varid ! NetCDF ID of a variable
     643integer :: ierr ! NetCDF routines return code
     644integer,dimension(4) :: id ! NetCDF IDs of the dimensions of the variable
     645integer,dimension(4) :: edges,corners
     646
     647#ifdef CPP_PARA
     648! Added to work in parallel mode
     649real dx3_glop(klon_glo,nsoilmx)
     650real dx3_glo(nbp_lon,nbp_lat,nsoilmx) ! to store a global 3D data set
     651real dx2_glop(klon_glo)
     652real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
     653real px2(ngrid)
     654#endif
     655
     656! 0. Do we ouput a diagsoil.nc file? If not just bail out now.
     657
     658! additional check: one can only output diagsoil.nc files
     659! in lon-lat case (or 1D)
     660if (grid_type==unstructured) then
     661  write(*,*) "writediagsoil: Error !!!"
     662  write(*,*) "diagsoil.nc outputs not possible on unstructured grids!!"
     663  call abort_physic("writediagsoil","impossible on unstructured grid",1)
     664endif
     665
     666! 1. Initialization step
     667if (firstname.eq."1234567890") then
     668  ! Store 'name' as 'firstname'
     669  firstname=name
     670  ! From now on, if 'name'.eq.'firstname', then it is a new time cycle
     671
     672  ! just to be sure, check that firstnom is large enough to hold nom
     673  if (len_trim(firstname).lt.len_trim(name)) then
     674    write(*,*) "writediagsoil: Error !!!"
     675    write(*,*) "   firstname string not long enough!!"
     676    write(*,*) "   increase its size to at least ",len_trim(name)
     677    call abort_physic("writediagsoil","firstname too short",1)
     678  endif
     679 
     680  ! Set output sample rate
     681  isample=int(ecritphy) ! same as for diagfi outputs
     682  ! Note ecritphy is known from control.h
     683 
     684  ! Create output NetCDF file
     685  if (is_master) then
     686   ierr=NF_CREATE(filename,NF_CLOBBER,nid)
     687   if (ierr.ne.NF_NOERR) then
     688    write(*,*)'writediagsoil: Error, failed creating file '//trim(filename)
     689    call abort_physic("writediagsoil","failed creating"//trim(filename),1)
     690   endif
     691  endif
     692
     693#ifdef CPP_PARA
     694   ! Gather inertiedat() soil thermal inertia on physics grid
     695   call Gather(inertiedat,inertiafi_glo)
     696   ! Gather cell_area() mesh area on physics grid
     697   call Gather(cell_area,areafi_glo)
     698#else
     699         inertiafi_glo(:,:)=inertiedat(:,:)
     700         areafi_glo(:)=cell_area(:)
     701#endif
     702
     703  if (is_master) then
     704   ! build inertia() and area()
     705   if (klon_glo>1) then
     706    do i=1,nbp_lon+1 ! poles
     707     inertia(i,1,1:nsoilmx)=inertiafi_glo(1,1:nsoilmx)
     708     inertia(i,nbp_lat,1:nsoilmx)=inertiafi_glo(klon_glo,1:nsoilmx)
     709     ! for area, divide at the poles by nbp_lon
     710     area(i,1)=areafi_glo(1)/nbp_lon
     711     area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
     712    enddo
     713    do j=2,nbp_lat-1
     714     ig0= 1+(j-2)*nbp_lon
     715     do i=1,nbp_lon
     716        inertia(i,j,1:nsoilmx)=inertiafi_glo(ig0+i,1:nsoilmx)
     717        area(i,j)=areafi_glo(ig0+i)
     718     enddo
     719     ! handle redundant point in longitude
     720     inertia(nbp_lon+1,j,1:nsoilmx)=inertia(1,j,1:nsoilmx)
     721     area(nbp_lon+1,j)=area(1,j)
     722    enddo
     723   endif
     724   
     725   ! write "header" of file (longitudes, latitudes, geopotential, ...)
     726   if (klon_glo>1) then ! general 3D case
     727     call iniwritesoil(nid,ngrid,inertia,area,nbp_lon+1,nbp_lat)
     728   else ! 1D model
     729     call iniwritesoil(nid,ngrid,inertiafi_glo(1,:),areafi_glo(1),1,1)
     730   endif
     731
     732  endif ! of if (is_master)
     733 
     734  ! set zitau to -1 to be compatible with zitau incrementation step below
     735  zitau=-1
     736 
     737else
     738  ! If not an initialization call, simply open the NetCDF file
     739  if (is_master) then
     740   ierr=NF_OPEN(filename,NF_WRITE,nid)
     741  endif
     742endif ! of if (firstname.eq."1234567890")
     743
     744! 2. Increment local time counter, if necessary
     745if (name.eq.firstname) then
     746  ! if we run across 'firstname', then it is a new time step
     747  zitau=zitau+iphysiq
     748  ! Note iphysiq is known from control.h
     749endif
     750
     751! 3. Write data, if the time index matches the sample rate
     752if (mod(zitau+1,isample).eq.0) then
     753
     754! 3.1 If first call at this date, update 'time' variable
     755  if (name.eq.firstname) then
     756    ntime=ntime+1
     757    date=float(zitau+1)/float(day_step)
     758    ! Note: day_step is known from control.h
     759   
     760    if (is_master) then
     761     ! Get NetCDF ID for "time"
     762     ierr=NF_INQ_VARID(nid,"time",varid)
     763     ! Add the current value of date to the "time" array
     764!#ifdef NC_DOUBLE
     765!     ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
     766!#else
     767     ierr=NF_PUT_VARA_REAL(nid,varid,[ntime],[1],[date])
     768!#endif
     769     if (ierr.ne.NF_NOERR) then
     770      write(*,*)"writediagsoil: Failed writing date to time variable"
     771      call abort_physic("writediagsoil","failed writing time",1)
     772     endif
     773    endif ! of if (is_master)
     774  endif ! of if (name.eq.firstname)
     775
     776! 3.2 Write the variable to the NetCDF file
     777if (dimpx.eq.3) then ! Case of a 3D variable
     778  ! A. Recast data along 'dynamics' grid
     779#ifdef CPP_PARA
     780  ! gather field on a "global" (without redundant longitude) array
     781  call Gather(px,dx3_glop)
     782!$OMP MASTER
     783  if (is_mpi_root) then
     784    call Grid1Dto2D_glo(dx3_glop,dx3_glo)
     785    ! copy dx3_glo() to dx3(:) and add redundant longitude
     786    data3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
     787    data3(nbp_lon+1,:,:)=data3(1,:,:)
     788  endif
     789!$OMP END MASTER
     790!$OMP BARRIER
     791#else
     792  if (klon_glo>1) then ! General case
     793   do l=1,nsoilmx
     794    ! handle the poles
     795    do i=1,nbp_lon+1
     796      data3(i,1,l)=px(1,l)
     797      data3(i,nbp_lat,l)=px(ngrid,l)
     798    enddo
     799    ! rest of the grid
     800    do j=2,nbp_lat-1
     801      ig0=1+(j-2)*nbp_lon
     802      do i=1,nbp_lon
     803        data3(i,j,l)=px(ig0+i,l)
     804      enddo
     805      data3(nbp_lon+1,j,l)=data3(1,j,l) ! extra (modulo) longitude
     806    enddo
     807   enddo
     808  else ! 1D model case
     809   data3_1d(1,1:nsoilmx)=px(1,1:nsoilmx)
     810  endif
     811#endif
     812 
     813  ! B. Write (append) the variable to the NetCDF file
     814  if (is_master) then
     815  ! B.1. Get the ID of the variable
     816  ierr=NF_INQ_VARID(nid,name,varid)
     817  if (ierr.ne.NF_NOERR) then
     818    ! If we failed geting the variable's ID, we assume it is because
     819    ! the variable doesn't exist yet and must be created.
     820    ! Start by obtaining corresponding dimensions IDs
     821    ierr=NF_INQ_DIMID(nid,"longitude",id(1))
     822    ierr=NF_INQ_DIMID(nid,"latitude",id(2))
     823    ierr=NF_INQ_DIMID(nid,"depth",id(3))
     824    ierr=NF_INQ_DIMID(nid,"time",id(4))
     825    ! Tell the world about it
     826    write(*,*) "====================="
     827    write(*,*) "writediagsoil: creating variable "//trim(name)
     828    call def_var(nid,name,title,units,4,id,varid,ierr)
     829  endif ! of if (ierr.ne.NF_NOERR)
     830 
     831  ! B.2. Prepare things to be able to write/append the variable
     832  corners(1)=1
     833  corners(2)=1
     834  corners(3)=1
     835  corners(4)=ntime
     836 
     837  if (klon_glo==1) then
     838    edges(1)=1
     839  else
     840    edges(1)=nbp_lon+1
     841  endif
     842  edges(2)=nbp_lat
     843  edges(3)=nsoilmx
     844  edges(4)=1
     845 
     846  ! B.3. Write the slab of data
     847!#ifdef NC_DOUBLE
     848!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3)
     849!#else
     850  if (klon_glo>1) then
     851    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3)
     852  else
     853    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3_1d)
     854  endif
     855!#endif
     856  if (ierr.ne.NF_NOERR) then
     857    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
     858               " to file "//trim(filename)//" at time",date
     859  endif
     860  endif ! of if (is_master)
     861
     862elseif (dimpx.eq.2) then ! Case of a 2D variable
     863
     864  ! A. Recast data along 'dynamics' grid
     865#ifdef CPP_PARA
     866  ! gather field on a "global" (without redundant longitude) array
     867  px2(:)=px(:,1)
     868  call Gather(px2,dx2_glop)
     869!$OMP MASTER
     870  if (is_mpi_root) then
     871    call Grid1Dto2D_glo(dx2_glop,dx2_glo)
     872    ! copy dx3_glo() to dx3(:) and add redundant longitude
     873    data2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
     874    data2(nbp_lon+1,:)=data2(1,:)
     875  endif
     876!$OMP END MASTER
     877!$OMP BARRIER
     878#else
     879  if (klon_glo>1) then ! general case
     880    ! handle the poles
     881    do i=1,nbp_lon+1
     882      data2(i,1)=px(1,1)
     883      data2(i,nbp_lat)=px(ngrid,1)
     884    enddo
     885    ! rest of the grid
     886    do j=2,nbp_lat-1
     887      ig0=1+(j-2)*nbp_lon
     888      do i=1,nbp_lon
     889        data2(i,j)=px(ig0+i,1)
     890      enddo
     891      data2(nbp_lon+1,j)=data2(1,j) ! extra (modulo) longitude
     892    enddo
     893  else ! 1D model case
     894    data2_1d=px(1,1)
     895  endif
     896#endif
     897
     898  ! B. Write (append) the variable to the NetCDF file
     899  if (is_master) then
     900  ! B.1. Get the ID of the variable
     901  ierr=NF_INQ_VARID(nid,name,varid)
     902  if (ierr.ne.NF_NOERR) then
     903    ! If we failed geting the variable's ID, we assume it is because
     904    ! the variable doesn't exist yet and must be created.
     905    ! Start by obtaining corresponding dimensions IDs
     906    ierr=NF_INQ_DIMID(nid,"longitude",id(1))
     907    ierr=NF_INQ_DIMID(nid,"latitude",id(2))
     908    ierr=NF_INQ_DIMID(nid,"time",id(3))
     909    ! Tell the world about it
     910    write(*,*) "====================="
     911    write(*,*) "writediagsoil: creating variable "//trim(name)
     912    call def_var(nid,name,title,units,3,id,varid,ierr)
     913  endif ! of if (ierr.ne.NF_NOERR)
     914
     915  ! B.2. Prepare things to be able to write/append the variable
     916  corners(1)=1
     917  corners(2)=1
     918  corners(3)=ntime
     919 
     920  if (klon_glo==1) then
     921    edges(1)=1
     922  else
     923    edges(1)=nbp_lon+1
     924  endif
     925  edges(2)=nbp_lat
     926  edges(3)=1
     927 
     928  ! B.3. Write the slab of data
     929!#ifdef NC_DOUBLE
     930!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2)
     931!#else
     932  if (klon_glo>1) then ! General case
     933    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2)
     934  else
     935    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,[data2_1d])
     936  endif
     937!#endif
     938  if (ierr.ne.NF_NOERR) then
     939    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
     940               " to file "//trim(filename)//" at time",date
     941  endif
     942  endif ! of if (is_master)
     943
     944elseif (dimpx.eq.0) then ! Case of a 0D variable
     945#ifdef CPP_PARA
     946  write(*,*) "writediagsoil: dimps==0 case not implemented in // mode!!"
     947  call abort_physic("writediagsoil","dimps==0 not implemented",1)
     948#endif
     949  ! A. Copy data value
     950  data0=px(1,1)
     951
     952  ! B. Write (append) the variable to the NetCDF file
     953  ! B.1. Get the ID of the variable
     954  ierr=NF_INQ_VARID(nid,name,varid)
     955  if (ierr.ne.NF_NOERR) then
     956    ! If we failed geting the variable's ID, we assume it is because
     957    ! the variable doesn't exist yet and must be created.
     958    ! Start by obtaining corresponding dimensions IDs
     959    ierr=NF_INQ_DIMID(nid,"time",id(1))
     960    ! Tell the world about it
     961    write(*,*) "====================="
     962    write(*,*) "writediagsoil: creating variable "//trim(name)
     963    call def_var(nid,name,title,units,1,id,varid,ierr)
     964  endif ! of if (ierr.ne.NF_NOERR)
     965
     966  ! B.2. Prepare things to be able to write/append the variable
     967  corners(1)=ntime
     968 
     969  edges(1)=1
     970
     971  ! B.3. Write the data
     972!#ifdef NC_DOUBLE
     973!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data0)
     974!#else
     975  ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,[data0])
     976!#endif
     977  if (ierr.ne.NF_NOERR) then
     978    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
     979               " to file "//trim(filename)//" at time",date
     980  endif
     981
     982endif ! of if (dimpx.eq.3) elseif (dimpx.eq.2) ...
     983endif ! of if (mod(zitau+1,isample).eq.0)
     984
     985! 4. Close the NetCDF file
     986if (is_master) then
     987  ierr=NF_CLOSE(nid)
     988endif
     989
     990end subroutine writediagsoilpem
     991
    581992
    582993END MODULE writediagpem_mod
Note: See TracChangeset for help on using the changeset viewer.