diff --git a/CaMa/src/cmf_calc_diag_mod.F90 b/CaMa/src/cmf_calc_diag_mod.F90 index 6648a75a..ad2fcc46 100644 --- a/CaMa/src/cmf_calc_diag_mod.F90 +++ b/CaMa/src/cmf_calc_diag_mod.F90 @@ -126,10 +126,10 @@ SUBROUTINE CMF_DIAG_AVERAGE D2RUNOFF_AVG(:,:) = D2RUNOFF_AVG(:,:) / DBLE(NADD) D2ROFSUB_AVG(:,:) = D2ROFSUB_AVG(:,:) / DBLE(NADD) IF ( LDAMOUT ) THEN - D2DAMINF_AVG(:,:) = D2DAMINF_AVG(:,:) / DBLE(NADD) + D2DAMINF_AVG(:,:) = D2DAMINF_AVG(:,:) / DBLE(NADD) ENDIF IF ( LWEVAP ) THEN - D2WEVAPEX_AVG(:,:) = D2WEVAPEX_AVG(:,:) / DBLE(NADD) + D2WEVAPEX_AVG(:,:) = D2WEVAPEX_AVG(:,:) / DBLE(NADD) ENDIF D1PTHFLW_AVG(:,:) = D1PTHFLW_AVG(:,:) /DBLE(NADD) END SUBROUTINE CMF_DIAG_AVERAGE diff --git a/CaMa/src/cmf_calc_outflw_mod.F90 b/CaMa/src/cmf_calc_outflw_mod.F90 index c2ad8013..b2986a27 100644 --- a/CaMa/src/cmf_calc_outflw_mod.F90 +++ b/CaMa/src/cmf_calc_outflw_mod.F90 @@ -267,9 +267,10 @@ SUBROUTINE CMF_CALC_INFLOW ENDDO !$OMP END PARALLEL DO - D2RIVINF(:,:)=P2RIVINF(:,:) !! needed for SinglePrecisionMode - D2FLDINF(:,:)=P2FLDINF(:,:) - + !D2RIVINF(:,:)=P2RIVINF(:,:) !! needed for SinglePrecisionMode + !D2FLDINF(:,:)=P2FLDINF(:,:) + D2RIVINF(1:NSEQALL,1)=P2RIVINF(1:NSEQALL,1) + D2FLDINF(1:NSEQALL,1)=P2FLDINF(1:NSEQALL,1) END SUBROUTINE CMF_CALC_INFLOW !#################################################################### END MODULE CMF_CALC_OUTFLW_MOD diff --git a/CaMa/src/cmf_ctrl_forcing_mod.F90 b/CaMa/src/cmf_ctrl_forcing_mod.F90 index f1e42b2f..99d37064 100644 --- a/CaMa/src/cmf_ctrl_forcing_mod.F90 +++ b/CaMa/src/cmf_ctrl_forcing_mod.F90 @@ -368,7 +368,7 @@ SUBROUTINE CMF_INPMAT_INIT_CDF(LECMF2LAKEC) !* Find levels dimension CALL NCERROR( NF90_INQUIRE_VARIABLE(NCID,VARID,dimids=VDIMIDS),'getting levI dimensions ') CALL NCERROR( NF90_INQUIRE_DIMENSION(NCID,VDIMIDS(1),len=INPNI),'getting time len ') - write(LOGNAM,*) 'Alocating INP*I: NXIN, NYIN, INPNI =', NXIN, NYIN, INPNI + write(LOGNAM,*) 'Allocating INP*I: NXIN, NYIN, INPNI =', NXIN, NYIN, INPNI allocate( INPXI(NXIN,NYIN,INPNI),INPYI(NXIN,NYIN,INPNI),INPAI(NXIN,NYIN,INPNI) ) write(LOGNAM,*)'INIT_MAP: inpaI:',TRIM(CINPMAT) @@ -550,9 +550,11 @@ SUBROUTINE CMF_FORCING_GET_CDF(PBUFF) IRECINP=INT( (KMIN-ROFCDF%NSTART)*60_JPIM,JPIM ) / INT(DTIN,JPIM) + 1 !! (second from netcdf start time) / (input time step) !*** 2. read runoff - CALL NCERROR( NF90_GET_VAR(ROFCDF%NCID,ROFCDF%NVARID(1),PBUFF(:,:,1),(/1,1,IRECINP/),(/NXIN,NYIN,1/)),'READING RUNOFF 1 ' ) + CALL NCERROR( NF90_GET_VAR(ROFCDF%NCID,ROFCDF%NVARID(1),PBUFF(:,:,1),& + (/1,1,IRECINP/),(/NXIN,NYIN,1/)),'READING RUNOFF 1 ' ) IF ( ROFCDF%NVARID(2) .NE. -1 ) THEN - CALL NCERROR( NF90_GET_VAR(ROFCDF%NCID,ROFCDF%NVARID(2),PBUFF(:,:,2),(/1,1,IRECINP/),(/NXIN,NYIN,1/)),'READING RUNOFF 2' ) + CALL NCERROR( NF90_GET_VAR(ROFCDF%NCID,ROFCDF%NVARID(2),PBUFF(:,:,2) & + ,(/1,1,IRECINP/),(/NXIN,NYIN,1/)),'READING RUNOFF 2' ) ENDIF write(LOGNAM,*) "CMF::FORCING_GET_CDF: read runoff:",IYYYYMMDD,IHHMM,IRECINP diff --git a/CaMa/src/cmf_ctrl_levee_mod.F90 b/CaMa/src/cmf_ctrl_levee_mod.F90 index a79586b1..2232b993 100644 --- a/CaMa/src/cmf_ctrl_levee_mod.F90 +++ b/CaMa/src/cmf_ctrl_levee_mod.F90 @@ -294,7 +294,7 @@ SUBROUTINE CMF_LEVEE_FLDSTG !$OMP PARALLEL DO REDUCTION(+:P0GLBSTOPRE2,P0GLBSTONEW2,P0GLBRIVSTO,P0GLBFLDSTO,P0GLBLEVSTO,P0GLBFLDARE) DO ISEQ=1, NSEQALL - +! DSTOALL = P2RIVSTO(ISEQ,1) + P2FLDSTO(ISEQ,1) + P2LEVSTO(ISEQ,1) DWTHINC = D2GRAREA(ISEQ,1) * D2RIVLEN(ISEQ,1)**(-1.) * DFRCINC !! width of each layer [m] IF( DSTOALL > P2RIVSTOMAX(ISEQ,1) )THEN diff --git a/CaMa/src/cmf_ctrl_maps_mod.F90 b/CaMa/src/cmf_ctrl_maps_mod.F90 index 3901b400..f8041376 100644 --- a/CaMa/src/cmf_ctrl_maps_mod.F90 +++ b/CaMa/src/cmf_ctrl_maps_mod.F90 @@ -44,7 +44,7 @@ MODULE CMF_CTRL_MAPS_MOD character(LEN=256) :: CRIVCLINC !! river map netcdf character(LEN=256) :: CRIVPARNC !! river parameter netcdf (WIDTH,HEIGHT, Manning, ground wateer delay) character(LEN=256) :: CMEANSLNC !! mean sea level netCDF - character(LEN=256) :: CMPIREGNC !! MPI Region netCDF + character(LEN=256) :: CMPIREGNC !! MPI region map in netcdf NAMELIST/NMAP/ CNEXTXY, CGRAREA, CELEVTN, CNXTDST, CRIVLEN, CFLDHGT, & CRIVWTH, CRIVHGT, CRIVMAN, CPTHOUT, CGDWDLY, CMEANSL, & @@ -313,8 +313,11 @@ SUBROUTINE CALC_REGION !! evenly allocate pixels to mpi nodes (updated in v4. ! integer(KIND=JPIM),SAVE :: IX,IY integer(KIND=JPIM),SAVE :: IREGION +#ifdef UseMPI_CMF #ifdef UseCDF_CMF - integer(KIND=JPIM) :: NCID,VARID +INTEGER(KIND=JPIM) :: NCID +INTEGER(KIND=JPIM) :: VARID +#endif #endif !$OMP THREADPRIVATE (IX) !================================================ @@ -381,7 +384,7 @@ SUBROUTINE CALC_REGION !! evenly allocate pixels to mpi nodes (updated in v4. write(LOGNAM,*) 'CALC_REGION: REGIONALL= ', REGIONALL write(LOGNAM,*) 'CALC_REGION: NSEQMAX=' , NSEQMAX - + WRITE(LOGNAM,*) 'CALC_REGION: NSEQALL=' , NSEQALL END SUBROUTINE CALC_REGION !========================================================== !+ diff --git a/CaMa/src/cmf_ctrl_mpi_mod.F90 b/CaMa/src/cmf_ctrl_mpi_mod.F90 index d987724d..d5a72b3a 100644 --- a/CaMa/src/cmf_ctrl_mpi_mod.F90 +++ b/CaMa/src/cmf_ctrl_mpi_mod.F90 @@ -1,5 +1,5 @@ MODULE CMF_CTRL_MPI_MOD -!! contains nothing is UseMPI_CMF is not defined +!! contains nothing if UseMPI_CMF is not defined #ifdef UseMPI_CMF !========================================================== !* PURPOSE: modules related to MPI usage @@ -19,7 +19,11 @@ MODULE CMF_CTRL_MPI_MOD ! See the License for the specific language governing permissions and limitations under the License. !========================================================== !** shared variables in module +#ifdef IFS_CMF + USE MPL_MODULE +#else USE MPI +#endif USE PARKIND1, only: JPIM, JPRB, JPRM, JPRD USE YOS_CMF_INPUT, only: LOGNAM USE YOS_CMF_MAP, only: REGIONALL, REGIONTHIS, MPI_COMM_CAMA @@ -37,11 +41,23 @@ MODULE CMF_CTRL_MPI_MOD ! -- CMF_DRV_END : Finalize CaMa-Flood ! !#################################################################### +#ifdef IFS_CMF +SUBROUTINE CMF_MPI_INIT(ICOMM_CMF) +IMPLICIT NONE + +INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: ICOMM_CMF +#else SUBROUTINE CMF_MPI_INIT IMPLICIT NONE +#endif !================================================ !*** 0. MPI specific setting REGIONTHIS=1 +#ifdef IFS_CMF + MPI_COMM_CAMA=ICOMM_CMF + REGIONALL=MPL_NPROC(MPI_COMM_CAMA) + REGIONTHIS=MPL_MYRANK(MPI_COMM_CAMA) +#else CALL MPI_Init(ierr) MPI_COMM_CAMA=MPI_COMM_WORLD @@ -50,7 +66,7 @@ SUBROUTINE CMF_MPI_INIT REGIONALL =Nproc REGIONTHIS=Nid+1 - +#endif ! For BUGFIX: Check MPI / OpenMPI is working or not. ! Write to standard output (log file is not opened yet) !!!!!#ifdef _OPENMP @@ -79,17 +95,23 @@ END SUBROUTINE CMF_MPI_END !#################################################################### SUBROUTINE CMF_MPI_AllReduce_R2MAP(R2MAP) - USE YOS_CMF_INPUT, only: NX,NY + USE YOS_CMF_INPUT, only: RMIS,NX,NY IMPLICIT NONE !* input/output real(KIND=JPRM),intent(inout) :: R2MAP(NX,NY) !* local variable real(KIND=JPRM) :: R2TMP(NX,NY) !================================================ - ! gather to master node - R2TMP(:,:)=1.E30 - CALL MPI_AllReduce(R2MAP,R2TMP,NX*NY,MPI_REAL4,MPI_MIN,MPI_COMM_CAMA,ierr) - R2MAP(:,:)=R2TMP(:,:) + ! gather to master node +#ifdef IFS_CMF + CALL MPL_ALLREDUCE(R2MAP,CDOPER='MIN',KCOMM=MPI_COMM_CAMA,KERROR=ierr) +#else + R2TMP(:,:)=RMIS + CALL MPI_AllReduce(R2MAP,R2TMP,NX*NY,MPI_REAL4,MPI_MIN,MPI_COMM_CAMA,ierr) + R2MAP(:,:)=R2TMP(:,:) +#endif + + END SUBROUTINE CMF_MPI_AllReduce_R2MAP !#################################################################### @@ -98,6 +120,7 @@ END SUBROUTINE CMF_MPI_AllReduce_R2MAP !#################################################################### SUBROUTINE CMF_MPI_AllReduce_R1PTH(R1PTH) + USE YOS_CMF_INPUT, ONLY: RMIS USE YOS_CMF_MAP, only: NPTHOUT, NPTHLEV IMPLICIT NONE !* input/output @@ -105,10 +128,14 @@ SUBROUTINE CMF_MPI_AllReduce_R1PTH(R1PTH) !* local variable real(KIND=JPRM) :: R1PTMP(NPTHOUT,NPTHLEV) !================================================ - ! gather to master node - R1PTMP(:,:)=1.E30 - CALL MPI_AllReduce(R1PTH,R1PTMP,NPTHOUT*NPTHLEV,MPI_REAL4,MPI_MIN,MPI_COMM_CAMA,ierr) - R1PTH(:,:)=R1PTMP(:,:) + ! gather to master node + #ifdef IFS_CMF + CALL MPL_ALLREDUCE(R1PTH,CDOPER='MIN',KCOMM=MPI_COMM_CAMA,KERROR=ierr) + #else + R1PTMP(:,:)=RMIS + CALL MPI_AllReduce(R1PTH,R1PTMP,NPTHOUT*NPTHLEV,MPI_REAL4,MPI_MIN,MPI_COMM_CAMA,ierr) + R1PTH(:,:)=R1PTMP(:,:) + #endif END SUBROUTINE CMF_MPI_AllReduce_R1PTH !#################################################################### @@ -116,7 +143,7 @@ END SUBROUTINE CMF_MPI_AllReduce_R1PTH !#################################################################### SUBROUTINE CMF_MPI_AllReduce_D2MAP(D2MAP) ! only used in netCDF restart file. (cannot be compiled due to a bug in MacOS mpif90) - USE YOS_CMF_INPUT, only: NX,NY + USE YOS_CMF_INPUT, only: DMIS, NX,NY IMPLICIT NONE !* input/output real(KIND=JPRB),intent(inout) :: D2MAP(NX,NY) @@ -124,13 +151,17 @@ SUBROUTINE CMF_MPI_AllReduce_D2MAP(D2MAP) real(KIND=JPRB) :: D2TMP(NX,NY) !================================================ ! gather to master node - D2TMP(:,:)=1.E30 +#ifdef IFS_CMF + CALL MPL_ALLREDUCE(D2MAP,CDOPER='MIN',KCOMM=MPI_COMM_CAMA,KERROR=ierr) +#else + D2TMP(:,:)=DMIS #ifdef SinglePrec_CMF CALL MPI_AllReduce(D2MAP,D2TMP,NX*NY,MPI_REAL4,MPI_MIN,MPI_COMM_CAMA,ierr) #else CALL MPI_AllReduce(D2MAP,D2TMP,NX*NY,MPI_REAL8,MPI_MIN,MPI_COMM_CAMA,ierr) #endif D2MAP(:,:)=D2TMP(:,:) +#endif END SUBROUTINE CMF_MPI_AllReduce_D2MAP !#################################################################### @@ -139,7 +170,7 @@ END SUBROUTINE CMF_MPI_AllReduce_D2MAP !#################################################################### SUBROUTINE CMF_MPI_AllReduce_P2MAP(P2MAP) ! only used in netCDF restart file. (cannot be compiled due to a bug in MacOS mpif90) - USE YOS_CMF_INPUT, only: NX,NY + USE YOS_CMF_INPUT, ONLY: DMIS, NX,NY IMPLICIT NONE !* input/output real(KIND=JPRD),intent(inout) :: P2MAP(NX,NY) @@ -147,9 +178,13 @@ SUBROUTINE CMF_MPI_AllReduce_P2MAP(P2MAP) real(KIND=JPRD) :: P2TMP(NX,NY) !================================================ ! gather to master node - P2TMP(:,:)=1.E30 + #ifdef IFS_CMF + CALL MPL_ALLREDUCE(P2MAP,CDOPER='MIN',KCOMM=MPI_COMM_CAMA,KERROR=ierr) + #else + P2TMP(:,:)=DMIS CALL MPI_AllReduce(P2MAP,P2TMP,NX*NY,MPI_REAL8,MPI_MIN,MPI_COMM_CAMA,ierr) P2MAP(:,:)=P2TMP(:,:) + #endif END SUBROUTINE CMF_MPI_AllReduce_P2MAP !#################################################################### @@ -158,7 +193,9 @@ END SUBROUTINE CMF_MPI_AllReduce_P2MAP !#################################################################### SUBROUTINE CMF_MPI_AllReduce_D1PTH(D1PTH) + USE YOS_CMF_INPUT, ONLY: DMIS USE YOS_CMF_MAP, only: NPTHOUT, NPTHLEV, PTH_UPST, PTH_DOWN + IMPLICIT NONE !* input/output real(KIND=JPRB),intent(inout) :: D1PTH(NPTHOUT,NPTHLEV) @@ -167,25 +204,29 @@ SUBROUTINE CMF_MPI_AllReduce_D1PTH(D1PTH) integer(KIND=JPIM) :: IPTH !================================================ ! gather to master node - DO IPTH=1,NPTHOUT - IF (PTH_UPST(IPTH)<=0 .or. PTH_DOWN(IPTH)<=0 ) THEN - D1PTH(IPTH,:)=1.E20 - ENDIF - ENDDO - D1PTMP(:,:)=1.E30 + DO IPTH=1,NPTHOUT + IF (PTH_UPST(IPTH)<=0 .OR. PTH_DOWN(IPTH)<=0 ) THEN + D1PTH(IPTH,:)=DMIS + ENDIF + END DO +#ifdef IFS_CMF + CALL MPL_ALLREDUCE(D1PTH,CDOPER='MIN',KCOMM=MPI_COMM_CAMA,KERROR=ierr) +#else + D1PTMP(:,:)=DMIS #ifdef SinglePrec_CMF - !! CALL MPI_Reduce(D1PTH,D1PTMP,NPTHOUT*NPTHLEV,MPI_REAL4,MPI_MIN,0,MPI_COMM_CAMA,ierr) CALL MPI_AllReduce(D1PTH,D1PTMP,NPTHOUT*NPTHLEV,MPI_REAL4,MPI_MIN,MPI_COMM_CAMA,ierr) #else CALL MPI_AllReduce(D1PTH,D1PTMP,NPTHOUT*NPTHLEV,MPI_REAL8,MPI_MIN,MPI_COMM_CAMA,ierr) #endif D1PTH(:,:)=D1PTMP(:,:) +#endif END SUBROUTINE CMF_MPI_AllReduce_D1PTH !#################################################################### !#################################################################### SUBROUTINE CMF_MPI_AllReduce_P1PTH(P1PTH) + USE YOS_CMF_INPUT, ONLY: DMIS USE YOS_CMF_MAP, only: NPTHOUT, NPTHLEV, PTH_UPST, PTH_DOWN IMPLICIT NONE !* input/output @@ -195,14 +236,18 @@ SUBROUTINE CMF_MPI_AllReduce_P1PTH(P1PTH) integer(KIND=JPIM) :: IPTH !================================================ ! gather to master node - DO IPTH=1,NPTHOUT - IF (PTH_UPST(IPTH)<=0 .or. PTH_DOWN(IPTH)<=0 ) THEN - P1PTH(IPTH,:)=1.E20 - ENDIF - ENDDO - P1PTMP(:,:)=1.E30 - CALL MPI_AllReduce(P1PTH,P1PTMP,NPTHOUT*NPTHLEV,MPI_REAL8,MPI_MIN,MPI_COMM_CAMA,ierr) - P1PTH(:,:)=P1PTMP(:,:) + DO IPTH=1,NPTHOUT + IF (PTH_UPST(IPTH)<=0 .OR. PTH_DOWN(IPTH)<=0 ) THEN + P1PTH(IPTH,:)=DMIS + ENDIF + END DO +#ifdef IFS_CMF + CALL MPL_ALLREDUCE(P1PTH,CDOPER='MIN',KCOMM=MPI_COMM_CAMA,KERROR=ierr) +#else + P1PTMP(:,:)=DMIS + CALL MPI_AllReduce(P1PTH,P1PTMP,NPTHOUT*NPTHLEV,MPI_REAL8,MPI_MIN,MPI_COMM_CAMA,ierr) + P1PTH(:,:)=P1PTMP(:,:) +#endif END SUBROUTINE CMF_MPI_AllReduce_P1PTH !#################################################################### @@ -219,8 +264,12 @@ SUBROUTINE CMF_MPI_ADPSTP(DT_MIN) !================================================ !*** MPI: use same DT in all node DT_LOC=DT_MIN - +#ifdef IFS_CMF + DT_OUT=DT_LOC + CALL MPL_ALLREDUCE(DT_OUT,CDOPER='MIN',KCOMM=MPI_COMM_CAMA,KERROR=ierr) +#else CALL MPI_AllReduce(DT_LOC, DT_OUT, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_CAMA,ierr) +#endif DT_MIN=DT_OUT write(LOGNAM,'(A,2F10.2)') "ADPSTP (MPI_AllReduce): DT_LOC->DTMIN", DT_LOC, DT_MIN diff --git a/CaMa/src/cmf_ctrl_nmlist_mod.F90 b/CaMa/src/cmf_ctrl_nmlist_mod.F90 index 103a4df9..77e703e1 100644 --- a/CaMa/src/cmf_ctrl_nmlist_mod.F90 +++ b/CaMa/src/cmf_ctrl_nmlist_mod.F90 @@ -226,8 +226,8 @@ SUBROUTINE CMF_CONFIG_NMLIST PMINSLP=1.E-5 !! minimum slope (kinematic wave) IMIS=-9999_JPIM - RMIS=-1.0E36_JPRM - DMIS=-1.0E36_JPRB + RMIS=1.E20_JPRM + DMIS=1.E20_JPRB CSUFBIN='.bin' CSUFVEC='.vec' diff --git a/CaMa/src/cmf_ctrl_vars_mod.F90 b/CaMa/src/cmf_ctrl_vars_mod.F90 index f3cfafaa..c06a8685 100644 --- a/CaMa/src/cmf_ctrl_vars_mod.F90 +++ b/CaMa/src/cmf_ctrl_vars_mod.F90 @@ -290,6 +290,7 @@ SUBROUTINE CMF_DIAG_INIT ALLOCATE(D2WINFILTEX_AVG(NSEQMAX,1)) D2WINFILTEX_AVG(:,:)=0._JPRB ENDIF + NADD=0 !*** 2b time-average 1D Diagnostics (bifurcation channel) diff --git a/CaMa/src/cmf_opt_outflw_mod.F90 b/CaMa/src/cmf_opt_outflw_mod.F90 index 09716162..adaf4655 100644 --- a/CaMa/src/cmf_opt_outflw_mod.F90 +++ b/CaMa/src/cmf_opt_outflw_mod.F90 @@ -292,7 +292,7 @@ SUBROUTINE CMF_CALC_OUTFLW_KINEMIX D2FLDOUT(ISEQ,1) = DARE_F * DVEL_F D2FLDOUT(ISEQ,1) = MIN( D2FLDOUT(ISEQ,1)*1._JPRD, P2FLDSTO(ISEQ,1)/DT ) ENDIF - ENDIF + ENDIF ENDDO !$OMP END PARALLEL DO @@ -371,7 +371,8 @@ SUBROUTINE CMF_CALC_OUTPRE DFLW = D2RIVDPH(ISEQ,1) DAREA = D2RIVWTH(ISEQ,1) * DFLW IF( DAREA>1.E-5 )THEN - D2RIVOUT_PRE(ISEQ,1) = DAREA * ( D2RIVMAN(ISEQ,1)**(-1.) * DFLW**(2./3.) * abs(DSLOPE)**(0.5) ) + D2RIVOUT_PRE(ISEQ,1) = DAREA * ( D2RIVMAN(ISEQ,1)**(-1.) * & + DFLW**(2./3.) * abs(DSLOPE)**(0.5) ) IF( DSLOPE<0._JPRB ) D2RIVOUT_PRE(ISEQ,1)=-D2RIVOUT_PRE(ISEQ,1) ELSE D2RIVOUT_PRE(ISEQ,1) = 0._JPRB @@ -381,7 +382,8 @@ SUBROUTINE CMF_CALC_OUTPRE DARE_F = P2FLDSTO(ISEQ,1) * D2RIVLEN(ISEQ,1)**(-1.) DARE_F = MAX( DARE_F - D2FLDDPH(ISEQ,1)*D2RIVWTH(ISEQ,1), 0._JPRB ) !! remove above river channel area IF( DARE_F>1.E-5 )THEN - D2FLDOUT_PRE(ISEQ,1) = DARE_F * ( PMANFLD**(-1.) * DFLW_F**(2./3.) * abs(DSLOPE_F)**(0.5) ) + D2FLDOUT_PRE(ISEQ,1) = DARE_F * ( PMANFLD**(-1.) * & + DFLW_F**(2./3.) * abs(DSLOPE_F)**(0.5)) IF( DSLOPE_F<0._JPRB ) D2FLDOUT_PRE(ISEQ,1)=-D2FLDOUT_PRE(ISEQ,1) ELSE D2FLDOUT_PRE(ISEQ,1) = 0._JPRB @@ -402,7 +404,8 @@ SUBROUTINE CMF_CALC_OUTPRE DFLW = MAX(DFLW,0._JPRB) IF( DFLW>1.E-5 )THEN - D1PTHFLW_PRE(IPTH,ILEV) = PTH_WTH(IPTH,ILEV) * DFLW * ( PTH_MAN(ILEV)**(-1.) * DFLW**(2./3.) * abs(DSLOPE)**(0.5) ) + D1PTHFLW_PRE(IPTH,ILEV) = PTH_WTH(IPTH,ILEV) * DFLW * & + ( PTH_MAN(ILEV)**(-1.) * DFLW**(2./3.) * abs(DSLOPE)**(0.5) ) IF( DSLOPE<0._JPRB ) D1PTHFLW_PRE(IPTH,ILEV)=-D1PTHFLW_PRE(IPTH,ILEV) ELSE D1PTHFLW_PRE(IPTH,ILEV) = 0._JPRB diff --git a/CaMa/src/cmf_utils_mod.F90 b/CaMa/src/cmf_utils_mod.F90 index e8e1b14d..96bb7c1e 100644 --- a/CaMa/src/cmf_utils_mod.F90 +++ b/CaMa/src/cmf_utils_mod.F90 @@ -534,9 +534,9 @@ END SUBROUTINE ENDIAN4I FUNCTION INQUIRE_FID() RESULT(FID) IMPLICIT NONE !* input/output - integer :: FID !< FILE ID + integer :: FID ! FILE ID !* local variable - logical :: I_OPENED !< FILE ID IS ALREADY USED or not? + logical :: I_OPENED ! FILE ID IS ALREADY USED or not? !================================================ DO FID = 10, 999 INQUIRE(FID,OPENED=I_OPENED) diff --git a/CaMa/src/parkind1.F90 b/CaMa/src/parkind1.F90 index e5c9c15d..27123d54 100644 --- a/CaMa/src/parkind1.F90 +++ b/CaMa/src/parkind1.F90 @@ -17,8 +17,8 @@ MODULE PARKIND1 !*** Integer Kinds integer,parameter :: JPIT = SELECTED_INT_KIND(2) integer,parameter :: JPIS = SELECTED_INT_KIND(4) - integer,parameter :: JPIM = SELECTED_INT_KIND(9) - integer,parameter :: JPIB = SELECTED_INT_KIND(12) + integer,parameter :: JPIM = SELECTED_INT_KIND(9) !! 4 byte integer + integer,parameter :: JPIB = SELECTED_INT_KIND(12) !! 8 byte long integer !Special integer type to be used for sensative adress calculations !should be *8 for a machine with 8byte adressing for optimum performance #ifdef ADDRESS64 @@ -30,15 +30,15 @@ MODULE PARKIND1 !*** Real Kinds integer,parameter :: JPRT = SELECTED_REAL_KIND(2,1) integer,parameter :: JPRS = SELECTED_REAL_KIND(4,2) - integer,parameter :: JPRM = SELECTED_REAL_KIND(6,37) +INTEGER, PARAMETER :: JPRM = SELECTED_REAL_KIND(6,37) !! 4 byte float #ifdef SinglePrec_CMF - integer,parameter :: JPRB = SELECTED_REAL_KIND(6,37) + integer,parameter :: JPRB = SELECTED_REAL_KIND(6,37) !! JPRB is switchable (4 byte in Single Precision Mode) #else integer,parameter :: JPRB = SELECTED_REAL_KIND(13,300) #endif ! Double real for C code and special places requiring ! higher precision. - integer,parameter :: JPRD = SELECTED_REAL_KIND(13,300) +INTEGER, PARAMETER :: JPRD = SELECTED_REAL_KIND(13,300) !! 8 byte double-precison float (primary used for precise water budget) !================================================ ! Logical Kinds for RTTOV.... diff --git a/CaMa/src/yos_cmf_diag.F90 b/CaMa/src/yos_cmf_diag.F90 index 52de8e88..05e12e43 100644 --- a/CaMa/src/yos_cmf_diag.F90 +++ b/CaMa/src/yos_cmf_diag.F90 @@ -46,7 +46,7 @@ MODULE YOS_CMF_DIAG !*** Average diagnostics real(KIND=JPRB),ALLOCATABLE,TARGET :: D2RIVOUT_AVG(:,:) !! average river discharge - real(KIND=JPRB),ALLOCATABLE,TARGET :: D2OUTFLW_AVG(:,:) !! average total outflow [m3/s] (rivout + fldout) !! bugfix v362 + real(KIND=JPRB),ALLOCATABLE,TARGET :: D2OUTFLW_AVG(:,:) !! average total outflow [m3/s] (rivout + fldout) real(KIND=JPRB),ALLOCATABLE,TARGET :: D2FLDOUT_AVG(:,:) !! average floodplain discharge real(KIND=JPRB),ALLOCATABLE,TARGET :: D2RIVVEL_AVG(:,:) !! average flow velocity real(KIND=JPRB),ALLOCATABLE,TARGET :: D2PTHOUT_AVG(:,:) !! flood pathway net outflow (2D) diff --git a/include/Makeoptions.gnu b/include/Makeoptions.gnu index ea970844..75b08623 100644 --- a/include/Makeoptions.gnu +++ b/include/Makeoptions.gnu @@ -44,7 +44,7 @@ CP = /bin/cp # DATM=-DNoAtom: activate when OMP ATOMIC calculation should be avoided (bit identical simulation) #---- #DMPI=-DUseMPI -DCDF=-DUseCDF +DCDF=-DUseCDF -DUseCDF_CMF #DATM=-DNoAtom CFLAGS=$(DMPI) $(DCDF) $(DATM) #---- diff --git a/main/MOD_RainSnowTemp.F90 b/main/MOD_RainSnowTemp.F90 index 2681d13d..32016378 100644 --- a/main/MOD_RainSnowTemp.F90 +++ b/main/MOD_RainSnowTemp.F90 @@ -109,11 +109,11 @@ SUBROUTINE rain_snow_temp (patchtype,& flfall = min(1.0_r8, max(0.0_r8,(forc_t - all_snow_t)*frac_rain_slope)) ELSEIF (trim(DEF_precip_phase_discrimination_scheme) == 'III') THEN - ! Phillip Harder and John Pomeroy (2013) - ! Estimating precipitation phase using a psychrometric energy - ! balance method . Hydrol Process, 27, 1901–1914 - ! Hydromet_Temp [K] - CALL Hydromet_Temp(forc_psrf,(forc_t-273.15),forc_q,t_hydro) + ! Phillip Harder and John Pomeroy (2013) + ! Estimating precipitation phase using a psychrometric energy + ! balance method . Hydrol Process, 27, 1901–1914 + ! Hydromet_Temp [K] + CALL hydromet_temp(forc_psrf,(forc_t-273.15),forc_q,t_hydro) IF(t_hydro > 3.0)THEN flfall = 1.0 ! fraction of liquid water within falling precip @@ -210,93 +210,98 @@ END SUBROUTINE NewSnowBulkDensity !!============================================== !----------------------------------------------------------------------------- - SUBROUTINE HYDROMET_TEMP(PPA, PTA, PQA,PTI) - !DESCRIPTION - !=========== - ! the temperature of a falling hydrometeor based on Harder, P., Pomeroy, J. (2013). - - !Original Author: - !------------------- + SUBROUTINE hydromet_temp(ppa, pta, pqa, pti) + ! DESCRIPTION + ! =========== + ! Computes the temperature of a falling hydrometeor based on Harder, P., Pomeroy, J. (2013). + + ! Original Author: + ! ---------------- ! V. Vionnet (11/2020) - !References: - !------------------- - !---Harder, P., Pomeroy, J. (2013). - ! Estimating precipitation phase using a psychrometric energy balance method - ! Hydrological Processes 27(13), 1901-1914. https://dx.doi.org/10.1002/hyp.9799 - !REVISION HISTORY - !---------------- - !---2023.07.30 Aobo Tan & Zhongwang Wei @ SYSU - - real(r8), intent(in) :: PPA ! Air pressure (Pa) - real(r8), intent(in) :: PTA ! Air temperature (deg C) - real(r8), intent(in) :: PQA ! Air specific humidity (kg/kg) - real(r8), intent(out) :: PTI ! Hydrometeo temprtature in deg C - real(r8) :: ZD ! diffusivity of water vapour in air [m^2 s-1] - real(r8) :: ZLAMBDAT ! thermal conductivity of air [J m^-1 s^-1 K^-1] - real(r8) :: ZL ! latent heat of sublimation of vaporisation[J kg^-1] - real(r8) :: ZRHODA ! density of dry air [kg m-3] - real(r8) :: ZRH ! relative humidity [-] - real(r8) :: RHO_VSAT_DIFF,ESAT,RHO_VSAT - real(r8) :: ZT,ZTINI,ZF,ZFDIFF,EVSAT - integer :: JITER - integer :: JJ,I,NN - - ! 1. Compute diffusivity of water vapour in air [m2 s-1] (Thorpe and Mason, 1966) - ZD = 2.063e-5 * ((PTA+273.15)/273.15)**1.75 - - ! 2. Compute thermal conductivity of air [J m-1 s-1 K-1] - ZLAMBDAT = 0.000063 * (PTA+273.15) + 0.00673 - - ! 3. Compute latent heat of sublimation or vaporisation (depending on air temperature) - IF(PTA <0.) THEN - ZL = 1000.0 * (2834.1 - 0.29 *PTA - 0.004*PTA**2.) - ELSE - ZL = 1000.0 * (2501.0 - (2.361 * PTA)) - ENDIF - - !TODO:check USE of dry air? - - ! 4. Compute density of dry air [kg m-3] - ZRHODA = PPA/(287.04*(PTA+273.15)) + + ! References: + ! ----------- + ! Harder, P., Pomeroy, J. (2013). + ! Estimating precipitation phase using a psychrometric energy balance method + ! Hydrological Processes 27(13), 1901-1914. https://dx.doi.org/10.1002/hyp.9799 + + ! REVISION HISTORY + ! ---------------- + ! 2023.07.30 Aobo Tan & Zhongwang Wei @ SYSU + + real(r8), intent(in) :: ppa ! Air pressure (Pa) + real(r8), intent(in) :: pta ! Air temperature (deg C) + real(r8), intent(in) :: pqa ! Air specific humidity (kg/kg) + real(r8), intent(out) :: pti ! Hydrometeor temperature in deg C + + real(r8) :: zd ! Diffusivity of water vapour in air [m^2 s^-1] + real(r8) :: zlambda ! Thermal conductivity of air [J m^-1 s^-1 K^-1] + real(r8) :: zl ! Latent heat of sublimation or vaporization [J kg^-1] + real(r8) :: zrhoda ! Density of dry air [kg m^-3] + real(r8) :: zrh ! Relative humidity [-] + real(r8) :: rho_vast_diff, esat, rho_vast + real(r8) :: zt, ztint, zf, zfdiff, evsat + integer :: JITER + integer :: JJ, I, NN + + ! 1. Compute diffusivity of water vapour in air [m^2 s^-1] (Thorpe and Mason, 1966) + zd = 2.063e-5 * ((pta + 273.15) / 273.15) ** 1.75 + + ! 2. Compute thermal conductivity of air [J m^-1 s^-1 K^-1] + zlambda = 0.000063 * (pta + 273.15) + 0.00673 + + ! 3. Compute latent heat of sublimation or vaporization (depending on air temperature) + IF (pta < 0.) THEN + zl = 1000.0 * (2834.1 - 0.29 * pta - 0.004 * pta ** 2.) + ELSE + zl = 1000.0 * (2501.0 - (2.361 * pta)) + END IF + + ! 4. Compute density of dry air [kg m^-3] + zrhoda = ppa / (287.04 * (pta + 273.15)) + ! 5. Compute saturated water vapour pressure [Pa] - IF(PTA>0) THEN - EVSAT = 611.0*EXP(17.27*PTA/(PTA+237.3)) + IF (pta > 0) THEN + evsat = 611.0 * EXP(17.27 * pta / (pta + 237.3)) ELSE - EVSAT = 611.0*EXP(21.87*PTA/(PTA+265.5)) - ENDIF - - ! 6. Solve iteratively to get Ti in Harder and Pomeroy (2013). using a Newton-Raphston approach - !set the 1st guess to PTA - ZT = PTA - !loop until convergence - DO JITER = 1,10 - ZTINI = ZT ! - - IF(ZT>0) THEN - ESAT = 611.0*EXP(17.27*ZT/(ZT+237.3)) + evsat = 611.0 * EXP(21.87 * pta / (pta + 265.5)) + END IF + + ! 6. Solve iteratively to get Ti in Harder and Pomeroy (2013) using a Newton-Raphson approach + ! Set the first guess to pta + zt = pta + + ! Loop until convergence + DO JITER = 1, 10 + ztint = zt + + IF (zt > 0) THEN + esat = 611.0 * EXP(17.27 * zt / (zt + 237.3)) ELSE - ESAT = 611.0*EXP(21.87*ZT/(ZT+265.5)) - ENDIF - - RHO_VSAT = ESAT/(461.5*(ZT+273.15)) ! Saturated water vapour density - - ZF = ZT - PTA - ZD*ZL/ZLAMBDAT * ( PQA*ZRHODA - RHO_VSAT) - - IF(ZT>0) THEN - RHO_VSAT_DIFF = 611.0/( 461.5*(ZT+273.15)) * EXP( 17.27*ZT/(ZT+ 237.3)) * & - (-1/(ZT+273.15) + 17.27* 237.3/((ZT+ 237.3))**2.) + esat = 611.0 * EXP(21.87 * zt / (zt + 265.5)) + END IF + + rho_vast = esat / (461.5 * (zt + 273.15)) ! Saturated water vapour density + + zf = zt - pta - zd * zl / zlambda * (pqa * zrhoda - rho_vast) + + IF (zt > 0) THEN + rho_vast_diff = 611.0 / (461.5 * (zt + 273.15)) * EXP(17.27 * zt / (zt + 237.3)) * & + (-1 / (zt + 273.15) + 17.27 * 237.3 / ((zt + 237.3) ** 2.)) ELSE - RHO_VSAT_DIFF = 611.0/( 461.5*(ZT+273.15)) * EXP( 21.87*ZT/(ZT+ 265.5)) * & - (-1/(ZT+273.15) + 21.87* 265.5/((ZT+ 265.5))**2.) - ENDIF - + rho_vast_diff = 611.0 / (461.5 * (zt + 273.15)) * EXP(21.87 * zt / (zt + 265.5)) * & + (-1 / (zt + 273.15) + 21.87 * 265.5 / ((zt + 265.5) ** 2.)) + END IF + + zfdiff = 1 + zd * zl / zlambda * rho_vast_diff + zt = ztint - zf / zfdiff + IF (ABS(zt - ztint) .LT. 0.01) EXIT + END DO + + pti = zt + + END SUBROUTINE hydromet_temp - ZFDIFF = 1 + ZD*ZL/ZLAMBDAT * RHO_VSAT_DIFF - ZT = ZTINI - ZF/ZFDIFF - IF(ABS(ZT- ZTINI) .lt. 0.01) EXIT - ENDDO - PTI = ZT - END SUBROUTINE HYDROMET_TEMP END MODULE MOD_RainSnowTemp