Develop2

README.md

@@ -32,7 +32,7 @@ Technical documentation freely available at http://dx.doi.org/10.5065/ew8g-yr95
 **Original HRLDAS model description paper**: Fei Chen, Kevin W. Manning, Margaret A. LeMone, Stanley B. Trier, Joseph G. Alfieri, Rita Roberts, Mukul Tewari, Dev Niyogi, Thomas W. Horst, Steven P. Oncley, Jeffrey B. Basara, and Peter D. Blanken, 2007: Description and Evaluation of the Characteristics of the NCAR High-Resolution Land Data Assimilation System. J. Appl. Meteor. Climatol., 46, 694–713.
 doi: http://dx.doi.org/10.1175/JAM2463.1
 
-**HRLDAS/Noah-MP version 5.0 model description paper**:  He, C., Valayamkunnath, P., Barlage, M., Chen, F., Gochis, D., Cabell, R., Schneider, T., Rasmussen, R., Niu, G.-Y., Yang, Z.-L., Niyogi, D., and Ek, M. (2023): Modernizing the open-source community Noah-MP land surface model (version 5.0) with enhanced modularity, interoperability, and applicability, Geosci. Model Dev., https://doi.org/10.5194/egusphere-2023-675, in review.
+**HRLDAS/Noah-MP version 5.0 model description paper**:  He, C., Valayamkunnath, P., Barlage, M., Chen, F., Gochis, D., Cabell, R., Schneider, T., Rasmussen, R., Niu, G.-Y., Yang, Z.-L., Niyogi, D., and Ek, M.: Modernizing the open-source community Noah with multi-parameterization options (Noah-MP) land surface model (version 5.0) with enhanced modularity, interoperability, and applicability, Geosci. Model Dev., 16, 5131–5151, https://doi.org/10.5194/gmd-16-5131-2023, 2023.
 
 
 ## HRLDAS GitHub structure

hrldas/HRLDAS_forcing/create_forcing.F

@@ -91,6 +91,7 @@ program create_forcing
   integer :: ihour
   integer :: numarg
   integer :: unmatched
+  logical :: ic_flag
 #ifndef _GFORTRAN_
   ! For some reason, gfortran does not like this declared external,
   ! but everyone else wants it declared external.  So we use #ifndef
@@ -441,8 +442,16 @@ program create_forcing
      ierr = nf90_close(ncid)
      call error_handler(ierr, "Problem closing Netcdf file")
 
+     ic_flag = .false.
      if ( full_ic_frq > -1 ) then  ! full_ic_frq == -1 turns off extra processing for initial conditions.
-        if ( ( (full_ic_frq>0) .and. (mod(ihour, full_ic_frq)==0) ) .or. ( ihour==0 ) ) then
+        if ( ihour == 0) ic_flag = .true. 
+        if ( full_ic_frq > 0 ) then
+           if (mod(ihour, full_ic_frq) == 0) then
+             ic_flag = .true.
+           endif
+        endif 
+
+        if ( ic_flag .eqv. .true. ) then
 
            call open_netcdf_for_output( trim(OutputDir)//"HRLDAS_setup_"//date(1:4)//date(6:7)//date(9:10)//date(12:13)//"_d"//hgrid, &
              ncid, version, geo_em%idim, geo_em%jdim, geo_em, .true.)
@@ -1259,7 +1268,11 @@ subroutine remap(datastruct, geo_em, buff)
           call latlon_to_ij(datastruct%proj, geo_em%lat(i,j), east_longitude, etax(i,j), etay(i,j))
 #else
           if (forcing_type=="ERA5") then
-            call latlon_to_ij(datastruct%proj, geo_em%lat(i,j), geo_em%lon(i,j)+360., etax(i,j), etay(i,j))
+            if (geo_em%lon(i,j) < 0) then
+              call latlon_to_ij(datastruct%proj, geo_em%lat(i,j), geo_em%lon(i,j)+360., etax(i,j), etay(i,j))
+            else
+              call latlon_to_ij(datastruct%proj, geo_em%lat(i,j), geo_em%lon(i,j), etax(i,j), etay(i,j))
+            endif
           else
             call latlon_to_ij(datastruct%proj, geo_em%lat(i,j), geo_em%lon(i,j), etax(i,j), etay(i,j))
           endif
@@ -2103,7 +2116,7 @@ subroutine read_era5_netcdf_unit(ncid, datastruct, ierr, current)
     array = [string('Tair_f_inst'),string('z'),string('lsm')] ! Only z (for geopotential) and lsm (for land-sea mask) are needed 
     nullify(datastruct%data)
     ! Match the netcdf attribute with one of the Vtable entries
-    VARSEARCH : do j = 0, vtable_count 
+    VARSEARCH : do j = 1, vtable_count 
     if (current%label == vtable(j)%name) then 
         VAR_NAME = array(j)%str
         datastruct%field  = vtable(j)%name
@@ -2198,7 +2211,7 @@ end function check_if_same_map
 !==============================================================================
 !==============================================================================
 
-end program create_forcing
+!end program create_forcing
 
 !==============================================================================
 !==============================================================================
@@ -2727,7 +2740,15 @@ subroutine interp_rainfall_nearest_neighbor(datastruct, newarr, mix, mjx, geo_em
 
         call latlon_to_ij(datastruct%proj, geo_em%lat(ii,jj), east_longitude, x, y)
 #else
-        call latlon_to_ij(datastruct%proj, geo_em%lat(ii,jj), geo_em%lon(ii,jj), x, y)
+       if (forcing_type=="ERA5") then
+          if (geo_em%lon(ii,jj) < 0) then
+             call latlon_to_ij(datastruct%proj, geo_em%lat(ii,jj), geo_em%lon(ii,jj)+360.0, x, y)
+          else
+             call latlon_to_ij(datastruct%proj, geo_em%lat(ii,jj), geo_em%lon(ii,jj), x, y)
+          endif
+       else
+          call latlon_to_ij(datastruct%proj, geo_em%lat(ii,jj), geo_em%lon(ii,jj), x, y)
+       endif
 #endif
 
         xn = nint(x)
@@ -3196,3 +3217,8 @@ end subroutine close_flnm
 !==============================================================================
 !==============================================================================
 
+end program create_forcing
+
+!==============================================================================
+!==============================================================================
+

hrldas/HRLDAS_forcing/lib/cio.c

@@ -69,7 +69,7 @@
 
 /* ****************************************************************** */
 
-copen(nunit, name, mode, err, oflag)
+int copen(nunit, name, mode, err, oflag)
  /*
   * nunit = UNIX file descriptor associated with file named *name*
   * name  = UNIX file name 
@@ -145,7 +145,7 @@ copen(nunit, name, mode, err, oflag)
 }
 
 /* ****************************************************************** */
-bnseek(fd, bread, mode, iprint)
+int bnseek(fd, bread, mode, iprint)
 
 /*  Move the read/write file pointer
        fd     : Unix file descriptor.
@@ -183,7 +183,7 @@ bnseek(fd, bread, mode, iprint)
 
 /* ****************************************************************** */
 
-bnread(fd, buf, nbuf, bread, ios, idiag)
+int bnread(fd, buf, nbuf, bread, ios, idiag)
  /*
   * fd = UNIX file descriptor number (NOT a Fortran unit) 
   * buf = area into which to read 
@@ -234,7 +234,7 @@ bnread(fd, buf, nbuf, bread, ios, idiag)
 
 /* ****************************************************************** */
 
-bnwrit(fd, buf, nbuf, bwritten, err, idiag)
+int bnwrit(fd, buf, nbuf, bwritten, err, idiag)
     int            *fd, *nbuf, buf[], *bwritten, *err, *idiag;
 
  /*
@@ -268,7 +268,7 @@ bnwrit(fd, buf, nbuf, bwritten, err, idiag)
 
 /* ****************************************************************** */
 
-cclose(nunit, iprint, err)
+int cclose(nunit, iprint, err)
 /*
 Close a C (UNIX?) file descriptor:
   nunit  : (INPUT)  : The C (UNIX?) file descriptor to close.

hrldas/HRLDAS_forcing/run/examples/ERA5/Prepare_ERA5-Land_Forcing.sh

@@ -0,0 +1,190 @@
+#!/bin/bash
+
+#################################################################################################
+#                                                                                               #
+#                             Alessandro Anav         JULY 2023                                 #
+#                                                                                               #
+#      This script extracts and prepares the variables needed to run NOAH-MP from ERA5-Land     # 
+#      The program expects one single file per year, namely data should span from               #
+#                  YYYY-01-01_00:00:00 to YYYY-12-31_23:00:00                                   #
+#                                                                                               #
+#                   For questions or issues: alessandro.anav@enea.it                            #
+#                                                                                               #
+#################################################################################################
+
+
+#################################################################################################
+#                                                                                               #
+#                                        USER DEFINED INPUT FILES                               # 
+#                                                                                               #
+#################################################################################################
+
+# set date for data processing
+inyear=2023
+inmonth=11
+inday=01
+preyear=2023
+premonth=10
+preday=31
+
+indate=${inyear}-${inmonth}-${inday}
+INFILE=../../../ERA5_forcing/ERA5-Land-Noah-MP-${indate}.grib
+predate=${preyear}-${premonth}-${preday}
+PREVIOUS_FILE=../../../ERA5_forcing/ERA5-Land-Noah-MP-${predate}.grib
+
+#################################################################################################
+#                                                                                               #
+#                                        INPUT VARIABLES                                        # 
+#                                                                                               #
+#################################################################################################
+
+vars_name=("10U" "10V" "2T" "2D" "SP" "SSRD" "STRD" "TP")
+vars_nout=("10U" "10V" "2T" "2D" "SP" "ACSSRD" "ACSTRD" "ACTP")
+vars_code=("var165" "var166" "var167" "var168" "var134" "var169" "var175" "var228")
+
+#################################################################################################
+#                                                                                               #
+#                                        START PROGRAM                                          # 
+#                                                                                               #
+#################################################################################################
+
+echo ""
+echo -e "\033[33m Processing Date:  ${indate}  FILE: $INFILE  \033[0m"
+echo ""
+
+#-------------------------------------------------------------------------- CREATING OUTPUT DIRECTORIES
+
+mkdir -p OUTPUTS
+
+for i in ${!vars_name[@]}; do
+  mkdir -p OUTPUTS/${vars_name[$i]}
+done
+
+mkdir -p OUTPUTS/E
+mkdir -p OUTPUTS/Q
+
+#################################################################################################
+#                                                                                               #
+#                                    PART 1: EXTRACT TIME-STEP VARIABLES                        # 
+#                                                                                               #
+#################################################################################################
+
+echo -e "\033[32m o-----------------> PART 1: EXTRACTING TIME-STEP VARIABLES \033[0m"
+
+#-------------------------------------------------------------------------- SPLIT THE VARIABLES FROM INPUT FILE
+grib_copy $INFILE   OUTPUTS/[shortName].grb
+
+for i in ${!vars_name[@]}; do
+   echo -e "\033[36m PROCESSING ${vars_nout[$i]} \033[0m"
+   var=${vars_name[$i]}
+   #-------------------------------------------------------------------------- SPLIT TIME
+   codes_split_file -1 OUTPUTS/${var,,}.grb
+
+   #-------------------------------------------------------------------------- LOOP OVER FILES TO RENAME THEM ACCORDING TO THE DATE 
+   for f1 in $(ls OUTPUTS/${var,,}.grb_*); do
+      date=(`cdo -s showtimestamp $f1`)
+      timestamp=$(sed -e 's/T/''/g' -e 's/:00:00/''/g' -e 's/-/''/g' <<< $date)
+      mv $f1 OUTPUTS/${vars_name[$i]}/ERA5_"${vars_nout[$i]}"_"${timestamp}"
+   done
+done
+echo ""
+
+#################################################################################################
+#                                                                                               #
+#                                    PART 2: DECUMULATE FIELDS                                  # 
+#                                                                                               #
+#################################################################################################
+
+echo -e "\033[32m o-----------------> PART 2: DECUMULATE FIELDS \033[0m"
+
+#-------------------------------------------------------------------------- EXTRACT THE ACCUMULATED VARIABLES FROM THE LAST TIME STEP OF PREVIOUS YEAR
+
+cdo -s -selvar,var169 -seldate,${predate}T23:00:00,${predate}T23:00:00 $PREVIOUS_FILE OUTPUTS/SSRD/ERA5_ACSSRD_"$preyear""$premonth""$preday"23
+cdo -s -selvar,var175 -seldate,${predate}T23:00:00,${predate}T23:00:00 $PREVIOUS_FILE OUTPUTS/STRD/ERA5_ACSTRD_"$preyear""$premonth""$preday"23
+cdo -s -selvar,var228 -seldate,${predate}T23:00:00,${predate}T23:00:00 $PREVIOUS_FILE OUTPUTS/TP/ERA5_ACTP_"$preyear""$premonth""$preday"23
+
+#-------------------------------------------------------------------------- DECUMULATE VARIABLES
+for i in ${!vars_nout[@]}; do
+   if [ ${vars_nout[$i]} = "ACSSRD" ] || [ ${vars_nout[$i]} = "ACSTRD" ] || [ ${vars_nout[$i]} = "ACTP" ] ; then
+
+      echo -e "\033[36m PROCESSING ${vars_nout[$i]} \033[0m"
+
+      FILES=OUTPUTS/${vars_name[$i]}/ERA5_"${vars_nout[$i]}"_"$inyear"*
+      prev_file="OUTPUTS/${vars_name[$i]}/ERA5_"${vars_nout[$i]}"_""$preyear""$premonth""$preday"23
+
+      for f in $FILES; do
+         # Find the hour reading the last two digits of the filename
+         file=$(basename "$f")
+         hour="${file: -2}"
+
+         # Create the output file name replacing the name of accumulated variable with instantaneous variable  
+         outfile=$(sed 's/'${vars_nout[$i]}'/'${vars_name[$i]}'/g' <<< $f)
+
+         # Decumulate the 02-24 fields and copy the 01 data from accumulated to instantaneous variable
+         if [[ "$hour" = "01" ]]; then
+            echo -e "\033[32m We are at 01: cp $f to $outfile  \033[0m"
+            cp $f $outfile
+         else
+            echo -e "\033[33m Decumulating: $f less $prev_file saved to $outfile  \033[0m"
+            cdo -s -sub $f $prev_file $outfile
+         fi
+         prev_file=$f
+      done
+      echo ""
+   fi
+done
+echo ""
+
+#################################################################################################
+#                                                                                               #
+#                               PART 3: COMPUTE SPECIFIC HUMIDITY                               # 
+#                                                                                               #
+#################################################################################################
+
+# Three steps are taken for calculating specific humidity
+# (1) calculate saturated water pressure e from Td 
+# (2) combine e and SP
+# (3) calculate q from e and SP
+
+echo -e "\033[32m o-----------------> PART 3: COMPUTING SPECIFIC HUMIDITY \033[0m"
+
+FILES=OUTPUTS/2D/ERA5_2D_"$inyear"*
+
+for f in $FILES; do
+   outfile1=$(sed 's/'2D'/'E'/g' <<< $f)
+   echo -e "\033[33m Computing saturated water pressure from $f and save results to $outfile1  \033[0m"
+   cdo -s -expr,"var1=6.112*exp(17.67*(var168-273.15)/(var168-273.15+243.5))" $f $outfile1
+
+   mkdir -p OUTPUTS/TMP
+   SP_file=$(sed 's/'2D'/'SP'/g' <<< $f)
+   outfile2=$(sed 's/'2D'/'TMP'/g' <<< $f)
+   echo -e "\033[33m Merging $outfile1 and $SP_file to  $outfile2 \033[0m"
+   cdo -s merge $outfile1 $SP_file $outfile2
+
+   outfile3=$(sed 's/'2D'/'Q'/g' <<< $f)
+   cdo -s -setparam,133.128 -expr,"var133=(0.622*var1)/(var134/100.-(0.378*var1))/1000." $outfile2 $outfile3
+done
+echo ""
+
+#################################################################################################
+#                                                                                               #
+#                                       CLEANING UP                                             # 
+#                                                                                               #
+#################################################################################################
+
+rm OUTPUTS/*.grb
+rm OUTPUTS/TP/ERA5_ACTP_"$inyear"*
+rm OUTPUTS/SSRD/ERA5_ACSSRD_"$inyear"*
+rm OUTPUTS/STRD/ERA5_ACSTRD_"$inyear"*
+rm -rf OUTPUTS/TMP
+rm -rf OUTPUTS/E
+
+#################################################################################################
+#                                                                                               #
+#                                        NICE ENDING MESSAGE                                    # 
+#                                                                                               #
+#################################################################################################
+
+echo -e "\033[32m o-----------------> ALL DONE  \033[0m"
+echo ""
+