Changeset 2975 for trunk/LMDZ.MARS/util/expandstartfi.F90

Timestamp:

May 30, 2023, 5:35:09 PM (21 months ago)

Author:

romain.vande

Message:

Mars util :

Adapt expandstartfi to subslope dimension + small correction regarding depreciated dimension "number_of_advected_field"

RV

File:

• trunk/LMDZ.MARS/util/expandstartfi.F90 (modified) (11 diffs)

trunk/LMDZ.MARS/util/expandstartfi.F90

@@ -20 +20 @@
 integer :: varid ! to store the ID of a variable
 integer :: datashape(4) ! to store dimension IDs of a given dataset
-integer :: corner(3),edges(3) ! to read data with a time axis
+integer :: corner(4),edges(4) ! to read data with a time axis
 character(len=90) :: varname ! name of a variable
 character(len=90) :: varatt ! name of attribute of a variable
@@ -28 +28 @@
 integer :: subsurface_layers_dimid
 integer :: nlayer_plus_1_dimid
+integer :: nslope_dimid
 integer :: number_of_advected_fields_dimid
 integer :: time_dimid
@@ -37 +38 @@
 integer :: subsurface_layers
 integer :: nlayer_plus_1
+integer :: nslope
 integer :: number_of_advected_fields
 integer :: timelen
 real,allocatable :: surf_field(:) ! to store a 1D field of physical_points elements
 real,allocatable :: subsurf_field(:,:) ! to store subsurface (2D field)
+real,allocatable :: subslope_field(:,:)
+real,allocatable :: nlayer_plus_1_field(:,:)
+real,allocatable :: subsurf_subslope_field(:,:,:)
 
 ! Output file dimensions:
@@ -51 +56 @@
 real,allocatable :: out_surf_field(:,:)
 real,allocatable :: out_subsurf_field(:,:,:)
+real,allocatable :: out_subslope_field(:,:,:)
+real,allocatable :: out_nlayer_plus_1_field(:,:,:)
+real,allocatable :: out_subsurf_subslope_field(:,:,:,:)
 ! Output file dimension IDs
 integer :: lon_dimid
 integer :: lat_dimid
 integer :: depth_dimid
+integer :: nslope_out_dimid
+integer :: nlayer_plus_1_out_dimid
 ! IDs of Output file variables
 integer :: lon_varid
@@ -143 +153 @@
 endif
 
-status=nf90_inq_dimid(inid,"number_of_advected_fields",number_of_advected_fields_dimid)
-if (status.ne.nf90_noerr) then
-  write(*,*)"Failed to find number_of_advected_fields dimension"
-  write(*,*)trim(nf90_strerror(status))
- stop
-endif
-status=nf90_inquire_dimension(inid,number_of_advected_fields_dimid,len=number_of_advected_fields)
-if (status.ne.nf90_noerr) then
-  write(*,*)"Failed to find number_of_advected_fields value"
+status=nf90_inq_dimid(inid,"nslope",nslope_dimid)
+if (status.ne.nf90_noerr) then
+  write(*,*)"Failed to find slope dimension, old startfi file"
+  nslope=0
+  write(*,*)trim(nf90_strerror(status))
+else
+status=nf90_inquire_dimension(inid,nslope_dimid,len=nslope)
+if (status.ne.nf90_noerr) then
+  write(*,*)"Failed to find nslope value"
   write(*,*)trim(nf90_strerror(status))
  stop
 else
-  write(*,*) " number_of_advected_fields = ",number_of_advected_fields
+  write(*,*) " nslope = ",nslope
+endif
 endif
 
@@ -177 +188 @@
 allocate(surf_field(physical_points))
 allocate(subsurf_field(physical_points,subsurface_layers))
+allocate(subslope_field(physical_points,nslope))
+allocate(nlayer_plus_1_field(physical_points,nlayer_plus_1))
+allocate(subsurf_subslope_field(physical_points,subsurface_layers,nslope))
 
 ! 2.1. Create output file
@@ -296 +310 @@
   stop
 endif
+! nslope:
+status=nf90_def_dim(outid,"nslope",nslope,nslope_out_dimid)
+if (status.ne.nf90_noerr) then
+  write(*,*) "Failed creating nslope dimension"
+  write(*,*)trim(nf90_strerror(status))
+  stop
+endif
+! nslope:
+status=nf90_def_dim(outid,"nlayer_plus_1",nlayer_plus_1,nlayer_plus_1_out_dimid)
+if (status.ne.nf90_noerr) then
+  write(*,*) "Failed creating nslope dimension"
+  write(*,*)trim(nf90_strerror(status))
+  stop
+endif
 
 !2.3. write variables to output file
@@ -379 +407 @@
 
 allocate(out_surf_field(lonlen,latlen))
-
-shape(:) = 0
+allocate(out_subsurf_field(lonlen,latlen,subsurface_layers))
+allocate(out_subslope_field(lonlen,latlen,nslope))
+allocate(out_nlayer_plus_1_field(lonlen,latlen,nlayer_plus_1))
+allocate(out_subsurf_subslope_field(lonlen,latlen,subsurface_layers,nslope))
+
 do ivar=1,nbinvars ! loop on all input variables
+  shape(:) = 0
   ! find out what dimensions are linked to this variable
   status=nf90_inquire_variable(inid,ivar,name=varname,ndims=nbdim,&
                                dimids=shape,natts=nbatt)
+
   if (((nbdim==1).and.(shape(1)==physical_points_dimid))&
   .or.((nbdim==2).and.(shape(1)==physical_points_dimid)&
@@ -473 +506 @@
     endif
   endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid))
-enddo ! of do i=1,nbinvars 
-
-! 2.4. 3D (subsurface) variables
-allocate(out_subsurf_field(lonlen,latlen,subsurface_layers))
-
-do ivar=1,nbinvars ! loop on all input variables
-  ! find out what dimensions are linked to this variable
-  status=nf90_inquire_variable(inid,ivar,name=varname,ndims=nbdim,&
-                               dimids=shape,natts=nbatt)
-  if (((nbdim==2).and.(shape(1)==physical_points_dimid)&
-                 .and.(shape(2)==subsurface_layers_dimid))&
-  .or.((nbdim==3).and.(shape(1)==physical_points_dimid)&
-                 .and.(shape(2)==subsurface_layers_dimid)&
-                 .and.(shape(3)==time_dimid))) then
+
+  if ((nbdim==2).and.(shape(1)==physical_points_dimid)&
+                 .and.(shape(2)==subsurface_layers_dimid)) then
     
     corner(1) = 1
@@ -570 +592 @@
     endif
   endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...)
+
+
+
+
+  if ((nbdim==3).and.(shape(1)==physical_points_dimid)&
+                 .and.(shape(2)==nslope_dimid)) then
+    
+    corner(1) = 1
+    corner(2) = 1
+    corner(3) = timelen
+    edges(1)  = physical_points
+    edges(2)  = nslope
+    edges(3)  = 1
+    
+    write(*,*) " processing: ",trim(varname)
+
+    ! load input data:
+    status=nf90_inq_varid(inid,varname,invarid)
+    status=nf90_get_var(inid,invarid,subslope_field,corner,edges)
+    
+    ! switch output file to to define mode
+    status=nf90_redef(outid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to swich to define mode"
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+    !define the variable
+    status=nf90_def_var(outid,trim(varname),NF90_REAL,&
+                         (/lon_dimid,lat_dimid,nslope_out_dimid/),varid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed creating variable ",trim(varname)
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+
+    ! variable attributes
+    do k=1,nbatt
+      status=nf90_inq_attname(inid,invarid,k,varatt)
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed getting attribute number",k," for ",trim(varname)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+      status=nf90_get_att(inid,invarid,trim(varatt),varattcontent)
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed loading attribute ",trim(varatt)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+
+      ! write the attribute and its value to output
+      status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent))
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed writing attribute ",trim(varatt)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+    enddo ! do k=1,nbatt
+
+    ! swich out of NetCDF define mode
+    status=nf90_enddef(outid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to swich out of define mode"
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+    
+    ! "convert" from subsurf_field(:,:) to out_subslope_field(:,:,:)
+    do i=1,lonlen
+      out_subslope_field(i,1,:)=subslope_field(1,:) ! North Pole
+      out_subslope_field(i,latlen,:)=subslope_field(physical_points,:) ! South Pole
+    enddo
+    do j=2,latlen-1
+      ig0=1+(j-2)*(lonlen-1)
+      do i=1,lonlen-1
+        out_subslope_field(i,j,:)=subslope_field(ig0+i,:)
+      enddo
+      out_subslope_field(lonlen,j,:)=out_subslope_field(1,j,:) ! redundant lon=180
+    enddo
+    
+    ! write the variable
+    status=nf90_put_var(outid,varid,out_subslope_field)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to write ",trim(varname)
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+  endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...)
+
+
+
+
+
+  if ((nbdim==3).and.(shape(1)==physical_points_dimid)&
+                 .and.(shape(2)==nlayer_plus_1_dimid)) then
+    
+    corner(1) = 1
+    corner(2) = 1
+    corner(3) = timelen
+    edges(1)  = physical_points
+    edges(2)  = nlayer_plus_1
+    edges(3)  = 1
+    
+    write(*,*) " processing: ",trim(varname)
+
+    ! load input data:
+    status=nf90_inq_varid(inid,varname,invarid)
+    status=nf90_get_var(inid,invarid,nlayer_plus_1_field,corner,edges)
+    
+    ! switch output file to to define mode
+    status=nf90_redef(outid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to swich to define mode"
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+    !define the variable
+    status=nf90_def_var(outid,trim(varname),NF90_REAL,&
+                         (/lon_dimid,lat_dimid,nlayer_plus_1_out_dimid/),varid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed creating variable ",trim(varname)
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+
+    ! variable attributes
+    do k=1,nbatt
+      status=nf90_inq_attname(inid,invarid,k,varatt)
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed getting attribute number",k," for ",trim(varname)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+      status=nf90_get_att(inid,invarid,trim(varatt),varattcontent)
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed loading attribute ",trim(varatt)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+
+      ! write the attribute and its value to output
+      status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent))
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed writing attribute ",trim(varatt)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+    enddo ! do k=1,nbatt
+
+    ! swich out of NetCDF define mode
+    status=nf90_enddef(outid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to swich out of define mode"
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+    
+    ! "convert" from subsurf_field(:,:) to out_nlayer_plus_1_field(:,:,:)
+    do i=1,lonlen
+      out_nlayer_plus_1_field(i,1,:)=nlayer_plus_1_field(1,:) ! North Pole
+      out_nlayer_plus_1_field(i,latlen,:)=nlayer_plus_1_field(physical_points,:) ! South Pole
+    enddo
+    do j=2,latlen-1
+      ig0=1+(j-2)*(lonlen-1)
+      do i=1,lonlen-1
+        out_nlayer_plus_1_field(i,j,:)=nlayer_plus_1_field(ig0+i,:)
+      enddo
+      out_nlayer_plus_1_field(lonlen,j,:)=out_nlayer_plus_1_field(1,j,:) ! redundant lon=180
+    enddo
+    
+    ! write the variable
+    status=nf90_put_var(outid,varid,out_nlayer_plus_1_field)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to write ",trim(varname)
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+  endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...)
+
+
+
+
+  if (nbdim==4) then
+    
+    corner(1) = 1
+    corner(2) = 1
+    corner(3) = 1
+    corner(4) = timelen
+    edges(1)  = physical_points
+    edges(2)  = subsurface_layers
+    edges(3)  = nslope
+    edges(4)  = 1
+    
+    write(*,*) " processing: ",trim(varname)
+
+    ! load input data:
+    status=nf90_inq_varid(inid,varname,invarid)
+    status=nf90_get_var(inid,invarid,subsurf_subslope_field,corner,edges)
+    
+    ! switch output file to to define mode
+    status=nf90_redef(outid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to swich to define mode"
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+    !define the variable
+    print *, "subsurface_layers_dimid", subsurface_layers_dimid
+    print *, "nslope_out_dimid", nslope_out_dimid
+    status=nf90_def_var(outid,trim(varname),NF90_REAL,&
+                         (/lon_dimid,lat_dimid,depth_dimid,nslope_out_dimid/),varid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed creating variable ",trim(varname)
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+
+    ! variable attributes
+    do k=1,nbatt
+      status=nf90_inq_attname(inid,invarid,k,varatt)
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed getting attribute number",k," for ",trim(varname)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+      status=nf90_get_att(inid,invarid,trim(varatt),varattcontent)
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed loading attribute ",trim(varatt)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+
+      ! write the attribute and its value to output
+      status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent))
+      if (status.ne.nf90_noerr) then
+        write(*,*) "Failed writing attribute ",trim(varatt)
+        write(*,*)trim(nf90_strerror(status))
+      stop
+      endif
+    enddo ! do k=1,nbatt
+
+    ! swich out of NetCDF define mode
+    status=nf90_enddef(outid)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to swich out of define mode"
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+    
+    ! "convert" from subsurf_field(:,:) to out_subsurf_subslope_field(:,:,:)
+    do i=1,lonlen
+      out_subsurf_subslope_field(i,1,:,:)=subsurf_subslope_field(1,:,:) ! North Pole
+      out_subsurf_subslope_field(i,latlen,:,:)=subsurf_subslope_field(physical_points,:,:) ! South Pole
+    enddo
+    do j=2,latlen-1
+      ig0=1+(j-2)*(lonlen-1)
+      do i=1,lonlen-1
+        out_subsurf_subslope_field(i,j,:,:)=subsurf_subslope_field(ig0+i,:,:)
+      enddo
+      out_subsurf_subslope_field(lonlen,j,:,:)=out_subsurf_subslope_field(1,j,:,:) ! redundant lon=180
+    enddo
+    
+    ! write the variable
+    status=nf90_put_var(outid,varid,out_subsurf_subslope_field)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "Failed to write ",trim(varname)
+      write(*,*)trim(nf90_strerror(status))
+      stop
+    endif
+  endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...)
+
+
+
 enddo ! of do i=1,nbinvars 
-
 
 ! 3. Finish things and cleanup
@@ -582 +877 @@
 endif
 
+  write(*,*)"Done writing file ",trim(outfile)
+
 end program expandstartfi