diff --git a/jobs/JRTOFS_GLO_NCODA_HYCOM_VAR b/jobs/JRTOFS_GLO_NCODA_3DVAR similarity index 100% rename from jobs/JRTOFS_GLO_NCODA_HYCOM_VAR rename to jobs/JRTOFS_GLO_NCODA_3DVAR diff --git a/parm/rtofs_glo.hycom.oanl.in b/parm/rtofs_glo.hycom.oanl.in index 596adde3..cd2f0d71 100755 --- a/parm/rtofs_glo.hycom.oanl.in +++ b/parm/rtofs_glo.hycom.oanl.in @@ -2,13 +2,9 @@ bias_lmt = 3., 8., bias_mld = .true., bias_opt = 'tsuv', - bias_wt = 1., - clm_wt = 0.5, 0.5, 1.0, 0.0 - cluster(2) = 2., - cluster(3) = 2., - cluster(5) = 2., - cluster(6) = 2., - conflict(1) = .false., + bias_wt = 0.5, + clm_wt = 0.5, 0.5, 1.2, 0.1, + cluster = 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, corr_mdl(2) = 1, debug(1) = .true., debug(2) = .true., @@ -28,17 +24,14 @@ fcst(3) = .true., fcst(4) = .true., fcst(6) = .true., + fcst_tau = 192, fgat = 1, 1, 1, 1, 1, 24, - fgat_rec = 8, 16, 17, 13, 14, 1, global = .true., - hafs = .false., himem = .true., ice_asm = .true., - ice_rec = 8, 15, 16, 1, lndz = 1.99, lvl_nmo = 2000., mask_opt = '3D', - mld_rec = 425, mld_src = 'modl', mx_asm_lvl = 2500., 2500., 2500., 2500., 2500., n_pass = 2, @@ -48,6 +41,7 @@ over(3) = 1., over(5) = 1., over(6) = 1., + priors = 'instant', prf_corr = 0.8, prf_hrs = 120, prf_slct(1) = 3., diff --git a/parm/rtofs_glo_ncoda_inc2mom6nc_lyr.input b/parm/rtofs_glo_ncoda_inc2mom6nc_lyr.input new file mode 100644 index 00000000..e7b01f33 --- /dev/null +++ b/parm/rtofs_glo_ncoda_inc2mom6nc_lyr.input @@ -0,0 +1,57 @@ +MOM.res.nc +MOM.res_1.nc +MOM.res_3.nc +MOM.res_4.nc +&TShincname +&UVincname +&IDM 'idm ' = longitudinal array size +&JDM 'jdm ' = latitudinal array size +&KDM 'kk ' = number of layers restart +&KDM 'kncoda' = number of NCODA levels equal to kk for layer increments + 0.50 'zl ' = nominal depth of layer 1 from MOM archive + 1.90 'zl ' = nominal depth of layer 2 + 4.42 'zl ' = nominal depth of layer 3 + 8.38 'zl ' = nominal depth of layer 4 + 13.185 'zl ' = nominal depth of layer 5 + 18.555 'zl ' = nominal depth of layer 6 + 24.895 'zl ' = nominal depth of layer 7 + 32.33 'zl ' = nominal depth of layer 8 + 40.33 'zl ' = nominal depth of layer 9 + 48.33 'zl ' = nominal depth of layer 10 + 56.33 'zl ' = nominal depth of layer 11 + 64.33 'zl ' = nominal depth of layer 12 + 72.33 'zl ' = nominal depth of layer 13 + 80.33 'zl ' = nominal depth of layer 14 + 88.33 'zl ' = nominal depth of layer 15 + 96.33 'zl ' = nominal depth of layer 16 + 104.33 'zl ' = nominal depth of layer 17 + 112.33 'zl ' = nominal depth of layer 18 + 120.33 'zl ' = nominal depth of layer 19 + 128.33 'zl ' = nominal depth of layer 20 + 136.33 'zl ' = nominal depth of layer 21 + 144.33 'zl ' = nominal depth of layer 22 + 152.33 'zl ' = nominal depth of layer 23 + 161.33 'zl ' = nominal depth of layer 24 + 174.53 'zl ' = nominal depth of layer 25 + 200.69 'zl ' = nominal depth of layer 26 + 239.84 'zl ' = nominal depth of layer 27 + 286.04 'zl ' = nominal depth of layer 28 + 340.56 'zl ' = nominal depth of layer 29 + 404.89 'zl ' = nominal depth of layer 30 + 480.8 'zl ' = nominal depth of layer 31 + 570.375 'zl ' = nominal depth of layer 32 + 676.075 'zl ' = nominal depth of layer 33 + 800.8 'zl ' = nominal depth of layer 34 + 947.975 'zl ' = nominal depth of layer 35 +1121.645 'zl ' = nominal depth of layer 36 +1326.57 'zl ' = nominal depth of layer 37 +1568.38 'zl ' = nominal depth of layer 38 +1899.27 'zl ' = nominal depth of layer 39 +2399.27 'zl ' = nominal depth of layer 40 +5449.635 'zl ' = nominal depth of layer 41 +&seatmpinc +&salintinc +&uvelinc +&vvelinc +NONE +&lyrthkname diff --git a/scripts/exrtofs_glo_ncoda_hycom_var.sh b/scripts/exrtofs_glo_ncoda_hycom_var.sh index b40ccb45..55a5a306 100755 --- a/scripts/exrtofs_glo_ncoda_hycom_var.sh +++ b/scripts/exrtofs_glo_ncoda_hycom_var.sh @@ -57,16 +57,18 @@ fi ln -sf $COMIN/ncoda/ocnqc $DATA # 1.b link in topo files -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.a ${DATA}/regional.grid.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.b ${DATA}/regional.grid.b -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.a ${DATA}/regional.depth.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.b ${DATA}/regional.depth.b -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.iso.sigma.a iso.sigma.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.iso.sigma.b iso.sigma.b - -# 1.c check if hycom restart file available -if [[ ! -s $COMINm1/${RUN}_${modID}.t00z.n00.restart.a || - ! -s $COMINm1/${RUN}_${modID}.t00z.n00.restart.b ]] +ln -f -s ${FIXrtofs}/regional.mom6.nc ${DATA}/regional.mom6.nc +ln -f -s ${FIXrtofs}/depth_GLBb0.08_09m11ob2_mom6.nc depth_GLBb0.08_09m11ob2_mom6.nc + +# 1.c check if restart files are available +if [[ ! -s $COMINm1/RESTART/${PDY}.000000.MOM.res.nc || + ! -s $COMINm1/RESTART/${PDY}.000000.MOM.res_1.nc || + ! -s $COMINm1/RESTART/${PDY}.000000.MOM.res_2.nc || + ! -s $COMINm1/RESTART/${PDY}.000000.MOM.res_3.nc || + ! -s $COMINm1/RESTART/${PDY}.000000.MOM.res_4.nc || + ! -s $COMINm1/RESTART/${PDY}.000000.MOM.res_5.nc || + ! -s $COMINm1/RESTART/${PDY}.000000.MOM.res_16.nc ]] + then $USHrtofs/${RUN}_abort.sh "FATAL ERROR: $job No restart file found" \ "in $COMINm1" 15 @@ -98,10 +100,9 @@ cat << eof2 > ogridnl kkm = 41, kko = 41, m = 4500, - n = 3298, + n = 3297, nnest = 1, nproj = -1, - rlat = 70.2, &end eof2 @@ -121,26 +122,26 @@ mkdir -p $log_dir echo timecheck RTOFS_GLO_HYCOM start setup at $(date) #NCODA setup -$EXECrtofs/rtofs_ncoda_setup 3D hycom ogridnl $ddtg > pout1 +$EXECrtofs/rtofs_ncoda_setup 3D mom ogridnl $ddtg > pout1 err=$?; export err ; err_chk echo " error from rtofs_ncoda_setup=",$err #NCODA prep echo timecheck RTOFS_GLO_HYCOM start prep at $(date) -mpiexec -n 72 --cpu-bind core $EXECrtofs/rtofs_ncoda_prep 3D hycom ogridnl $ddtg > pout2 +mpiexec -n 72 --cpu-bind core $EXECrtofs/rtofs_ncoda_prep 3D mom ogridnl $ddtg > pout2 err=$?; export err ; err_chk echo " error from rtofs_ncoda_prep=",$err #NCODA var echo timecheck RTOFS_GLO_HYCOM start ncoda3d at $(date) -mpiexec -n $NPROCS --cpu-bind core $EXECrtofs/rtofs_ncoda 3D hycom ogridnl $ddtg > pout3 +mpiexec -n $NPROCS --cpu-bind core $EXECrtofs/rtofs_ncoda 3D mom ogridnl $ddtg > pout3 err=$?; export err ; err_chk echo " error from rtofs_ncoda=",$err #NCODA post echo timecheck RTOFS_GLO_HYCOM start post at $(date) . prep_step -mpiexec -n $NPROCS --cpu-bind core $EXECrtofs/rtofs_ncoda_post 3D hycom ogridnl $ddtg relax > pout4 +mpiexec -n $NPROCS --cpu-bind core $EXECrtofs/rtofs_ncoda_post 3D mom ogridnl $ddtg relax > pout4 err=$?; export err ; err_chk echo " error from rtofs_ncoda_post=",$err @@ -148,18 +149,13 @@ echo " error from rtofs_ncoda_post=",$err [[ -f fort.32 ]] && mv fort.32 $log_dir/hycom_var.$ddtg.rej [[ -f fort.33 ]] && mv fort.33 $log_dir/hycom_var.$ddtg.prf [[ -f fort.34 ]] && mv fort.34 $log_dir/hycom_var.$ddtg.gpt -[[ -f fort.36 ]] && mv fort.36 $log_dir/hycom_var.$ddtg.mvo [[ -f fort.37 ]] && mv fort.37 $log_dir/hycom_var.$ddtg.drc -[[ -f fort.38 ]] && mv fort.38 $log_dir/hycom_var.$ddtg.lyp [[ -f fort.39 ]] && mv fort.39 $log_dir/hycom_var.$ddtg.fix -[[ -f fort.40 ]] && mv fort.40 $log_dir/hycom_var.$ddtg.sus [[ -f fort.41 ]] && mv fort.41 $log_dir/hycom_var.$ddtg.dup [[ -f fort.42 ]] && mv fort.42 $log_dir/hycom_var.$ddtg.ssh [[ -f fort.52 ]] && mv fort.52 $log_dir/hycom_var.$ddtg.sal [[ -f fort.67 ]] && mv fort.67 $log_dir/hycom_var.$ddtg.obs [[ -f fort.68 ]] && mv fort.68 $log_dir/hycom_var.$ddtg.grd -[[ -f fort.69 ]] && mv fort.69 $log_dir/hycom_var.$ddtg.via -[[ -f fort.88 ]] && mv fort.88 $log_dir/hycom_var.$ddtg.dbg # create data coverage graphics DoGraphics=NO diff --git a/scripts/exrtofs_glo_ncoda_inc.sh b/scripts/exrtofs_glo_ncoda_inc.sh index 70960c8c..2abe2240 100755 --- a/scripts/exrtofs_glo_ncoda_inc.sh +++ b/scripts/exrtofs_glo_ncoda_inc.sh @@ -15,6 +15,7 @@ set -xa # Script history log: # # 2020-07-30 Dan Iredell # # 2023-02-08 Dmitry Dukhovskoy modified for updated ncoda_archv_lyrinc # +# 2026-01-12 Zulema Garraffo modified for MOM # # ############################################################################### @@ -34,34 +35,21 @@ incup_hours=6 dtg=${PDYm1}00 dtgm1=$($EXECrtofs/rtofs_dtg $dtg -d -1) dtgm2=$($EXECrtofs/rtofs_dtg $dtg -h -$incup_hours) -hday=$($USHrtofs/rtofs_date_normal2hycom.sh $dtg) -hday2=$(echo $hday $incup_hours 24 | awk '{printf("%9.3f", $1-($2/$3))}') -jday=$($USHutil/date2jday.sh ${dtg:0:8}) -jday2=$($USHutil/date2jday.sh ${dtgm2:0:8}) -archday=${jday:0:4}_${jday:4:3}_${dtg:8:2} -archday2=${jday2:0:4}_${jday2:4:3}_${dtgm2:8:2} - +jday=$($EXECrtofs/rtofs_dtg -f %Y%j $dtg) +fcst='0024' echo dtg12 $dtg $dtgm1 $dtgm2 -echo hday $hday $hday2 jday $jday $jday2 -echo archday $archday $archday2 mode=incup -BLKDATA_FILE=${PARMrtofs}/${RUN}_${modID}.${inputgrid}.${mode}.blkdat.input -IDM=$(cat ${BLKDATA_FILE} | grep idm | cut -d' ' -f1 | tr -d '[:space:]') -JDM=$(cat ${BLKDATA_FILE} | grep jdm | cut -d' ' -f1 | tr -d '[:space:]') -JDMA=$(expr ${JDM} \- 1) +#inputgrid=0.08, change to 0p08 +reg=GLB +DEPTH_FILE=${FIXrtofs}/depth_${reg}${inputgrid}_09m11ob2_mom6.nc +IDM=$(ncdump -h ${DEPTH_FILE} | grep 'nx =' | cut -d' ' -f3) +JDM=$(ncdump -h ${DEPTH_FILE} | grep 'ny =' | cut -d' ' -f3) +KDM=41 # or get from restart file? SIZN="${IDM}x${JDM}" -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.a regional.grid.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.b regional.grid.b -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.a regional.depth.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.b regional.depth.b -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.iso.sigma.a iso.sigma.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.iso.sigma.b iso.sigma.b -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.tbaric.a tbaric.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.tbaric.b tbaric.b -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.relax_ssh.a relax.ssh.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.relax_ssh.b relax.ssh.b +ln -f -s ${FIXrtofs}/depth_${reg}.${inputgrid}_09m11ob2_mom6.nc . +ln -f -s ${FIXrtofs}/regional.mom6.nc . # 2. link to ncoda hycom var restart files @@ -69,23 +57,23 @@ typet=seatmp_lyr_1o${SIZN} types=salint_lyr_1o${SIZN} typeu=uucurr_lyr_1o${SIZN} typev=vvcurr_lyr_1o${SIZN} -typep=lyrprs_lyr_1o${SIZN} -typec=icecov_sfc_1o${SIZN} - -ln -sf $COMINm1/rtofs_glo.t00z.n00.archv.a archv.${archday}.a -ln -sf $COMINm1/rtofs_glo.t00z.n00.archv.b archv.${archday}.b +typethbg=lyrthk_lyr_1o${SIZN} -export lyrprinc=lyrprs_${dtg}_analinc export salininc=salint_${dtg}_analinc -export stempinc=seatmp_${dtg}_analinc -export upvelinc=upvel_${dtg}_analinc # u vel increm on p-grid -export vpvelinc=vpvel_${dtg}_analinc +export tempinc=seatmp_${dtg}_analinc +export uvelinc=uvel_${dtg}_analinc # u vel increm on p-grid +export vvelinc=vvel_${dtg}_analinc +export lyrthbg=lyrthk_${dtgm1}_fcstfld + +#names can be changed. Names in INPUT will be MOM.inc.TSzh.nc, MOM.inc.UV.nc +export TShar=MOM.res_Y${jday:0:4}_D${jday:4:3}_S00000_inc.TSzh.nc +export UVar=MOM.res_Y${jday:0:4}_D${jday:4:3}_S00000_inc.TSzh.nc # Check for the existence of analysis increment files # These are needed to create the HYCOM incremental update file # Temperature if [ -e $COMIN/ncoda/hycom_var/restart/${typet}_${dtg}_0000_analinc ]; then - ln -sf $COMIN/ncoda/hycom_var/restart/${typet}_${dtg}_0000_analinc ./${stempinc} + ln -sf $COMIN/ncoda/hycom_var/restart/${typet}_${dtg}_0000_analinc ./${tempinc} else msg="$COMIN/ncoda/hycom_var/restart/${typet}_${dtg}_0000_analinc is missing" err_exit $msg @@ -99,99 +87,74 @@ else fi # Current - U-component if [ -e $COMIN/ncoda/hycom_var/restart/${typeu}_${dtg}_0000_analinc ]; then - ln -sf $COMIN/ncoda/hycom_var/restart/${typeu}_${dtg}_0000_analinc ./${upvelinc} + ln -sf $COMIN/ncoda/hycom_var/restart/${typeu}_${dtg}_0000_analinc ./${uvelinc} else msg="$COMIN/ncoda/hycom_var/restart/${typeu}_${dtg}_0000_analinc is missing" err_exit $msg fi # Current - V-component if [ -e $COMIN/ncoda/hycom_var/restart/${typev}_${dtg}_0000_analinc ]; then - ln -sf $COMIN/ncoda/hycom_var/restart/${typev}_${dtg}_0000_analinc ./${vpvelinc} + ln -sf $COMIN/ncoda/hycom_var/restart/${typev}_${dtg}_0000_analinc ./${vvelinc} else msg="$COMIN/ncoda/hycom_var/restart/${typev}_${dtg}_0000_analinc is missing" err_exit $msg fi -# Layer Pressure -if [ -e $COMIN/ncoda/hycom_var/restart/${typep}_${dtg}_0000_analinc ]; then - ln -sf $COMIN/ncoda/hycom_var/restart/${typep}_${dtg}_0000_analinc ./${lyrprinc} +# Background state layer thickness +if [ -e $COMIN/ncoda/hycom_var/restart/${typethbg}_${dtgm1}_0024_fcstfld ]; then + ln -sf $COMIN/ncoda/hycom_var/restart/${typethbg}_${dtgm1}_0024_fcstfld ./${lyrthbg} else - msg="$COMIN/ncoda/hycom_var/restart/${typep}_${dtg}_0000_analinc is missing" + msg="$COMIN/ncoda/hycom_var/restart/${typethbg}_${dtgm1}_0024_fcstfld is missing" err_exit $msg fi -# Ice Coverage -if [ -e $COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld ]; then - ln -sf $COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld ./icecov_${dtg}_analfld + +# Link the MOM6 template restart files +if [ -e $COMINm1/RESTART/${dtg:0:8}_000000.MOM.res.nc]; then + ln -sf $COMINm1/RESTART/${dtg:0:8}_000000.MOM.res.nc MOM.res.nc else - msg="$COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld is missing" + msg="$COMINm1/${dtg:0:8}_000000.MOM.res.nc is missing" + err_exit $msg +fi +if [ -e $COMINm1/RESTART/${dtg:0:8}_000000.MOM.res_1.nc]; then + ln -sf $COMINm1/RESTART/${dtg:0:8}_000000.MOM_1.res.nc MOM.res.nc +else + msg="$COMINm1/${dtg:0:8}_000000.MOM_1.res.nc is missing" + err_exit $msg +fi +if [ -e $COMINm1/RESTART/${dtg:0:8}_000000.MOM_3.res.nc]; then + ln -sf $COMINm1/RESTART/${dtg:0:8}_000000.MOM_3.res.nc MOM.res.nc +else + msg="$COMINm1/${dtg:0:8}_000000.MOM_3.res.nc is missing" + err_exit $msg +fi +if [ -e $COMINm1/RESTART/${dtg:0:8}_000000.MOM_4.res.nc]; then + ln -sf $COMINm1/RESTART/${dtg:0:8}_000000.MOM_4.res.nc MOM.res.nc +else + msg="$COMINm1/${dtg:0:8}_000000.MOM_4.res.nc is missing" err_exit $msg fi -#create ssmi.r file -rm -f ssmi1.a ssmi2.a ssmi.$dtg.r -$EXECrtofs/rtofs_raw2hycom icecov_${dtg}_analfld $IDM $JDM 999.00 ssmi1.a > ssmi1.b -err=$?; export err ; err_chk -echo " error from rtofs_raw2hycom=",$err -$EXECrtofs/rtofs_hycom_expr ssmi1.a "ONE" $IDM $JDM 0.01 0 ssmi2.a > ssmi2.b -err=$?; export err ; err_chk -echo " error from rtofs_hycom_expr=",$err -$EXECrtofs/rtofs_hycom2raw8 ssmi2.a $IDM $JDM 1 1 $IDM $JDMA ssmi.$dtg.r -err=$?; export err ; err_chk -echo " error from rtofs_hycom2raw8=",$err - -ar=archv_1_inc.${archday} -rm -f $ar.[a,b] # copy modify input file with local vars -cp ${PARMrtofs}/${RUN}_${modID}.ncoda_archv_lyr.input ./ncoda_archv.input -sed -i -e "s/&archday/$archday/" \ - -e "s/&archname/$ar/" \ +cp ${PARMrtofs}/${RUN}_${modID}.ncoda_inc2mom6nc.input ./ncoda_inc2mom6nc.input +sed -i -e "s/&TShincname/$TShar/" \ + -e "s/&UVincname/$UVar/" \ -e "s/&IDM/$IDM/g" \ -e "s/&JDM/$JDM/g" \ - -e "s/&dtg/$dtg/g" \ - -e "s/&lyrprsinc/${lyrprinc}/g" \ + -e "s/&KDM/$KDM/g" \ + -e "s/&seatmpinc/${tempinc}/g" \ -e "s/&salintinc/${salininc}/g" \ - -e "s/&seatmpinc/${stempinc}/g" \ - -e "s/&upvelinc/${upvelinc}/g" \ - -e "s/&vpvelinc/${vpvelinc}/g" ./ncoda_archv.input -#ln -s ${PARMrtofs}/${RUN}_${modID}.zlevels zi.txt + -e "s/&uvelinc/${uvelinc}/g" \ + -e "s/&vvelinc/${vvelinc}/g" \ + -e "s/&lyrthknam/${lyrthbg}" \ -$EXECrtofs/rtofs_ncoda_archv_lyrinc < ncoda_archv.input >> $pgmout +$EXECrtofs/rtofs_ncodaz_inc2mom6nc_glb_lyr.x < ncoda_inc2mom6nc_lyr.input >> $pgmout err=$?; export err ; err_chk -echo " error from rtofs_ncoda_archv_lyrinc=",$err +echo " error from rtofs_ncodaz_inc2mom6nc_glb_lyr=",$err -# -# calculate increment file for assimilation -# -/bin/rm -f archd_1.${dtg}.a archd_1.${dtg}.b -cat << E-o-D > hycom_diff.input -41 'kk ' = number of layers involved -1.0 'nscale' = scale difference by 1.0/nscale -0 'nbox ' = smooth difference over a 2*nbox+1 square -archv_1_inc.${archday}.a -archv.${archday}.a -archd_1.${dtg} -Analysis - Background -17T Sigma2*; GDEM4.2; KPP; SeaWiFS chl; HYCOM+CICE; A=20;Smag=.05; -Z(7):1-7,Z(16):8,Z(2):10-16,Z(13):dp00/f/x=36/1.18/262;Z(3):400-600m; GPCPsnow -sigma:84-14m; depth_GLBb0.08_11; apply offlux to CICE; 2.2.99DHi-900 -E-o-D -date -$EXECrtofs/rtofs_hycom_diff < hycom_diff.input -err=$?; export err ; err_chk -echo " error from rtofs_hycom_diff=",$err +msg="THE RTOFS_GLO_NCODA_INC JOB HAS ENDED NORMALLY on $(hostname) at $(date)" +postmsg "$msg" # -# change time on archd*.b file +# calculate increment file for assimilation # -sed -e "s/${hday}/${hday2}/g" archd_1.${dtg}.b > archd_1.${dtg}.b2 - -cp archd_1.$dtg.a $COMOUT/rtofs_glo.incupd.$archday2.a -cp archd_1.$dtg.b2 $COMOUT/rtofs_glo.incupd.$archday2.b -cp archv_1_inc.$archday.a $COMOUT/rtofs_glo.archv_1_inc.$archday.a -cp archv_1_inc.$archday.b $COMOUT/rtofs_glo.archv_1_inc.$archday.b -cp ssmi.$dtg.r $COMOUT/rtofs_glo.ssmi.$dtg.r - -msg="THE RTOFS_GLO_NCODA_INC JOB HAS ENDED NORMALLY on $(hostname) at $(date)" -postmsg "$msg" - diff --git a/scripts/exrtofs_glo_ncoda_qc.sh b/scripts/exrtofs_glo_ncoda_qc.sh index 7bdc6a0f..c6d79b6a 100755 --- a/scripts/exrtofs_glo_ncoda_qc.sh +++ b/scripts/exrtofs_glo_ncoda_qc.sh @@ -83,11 +83,11 @@ for v in glbl_var hycom_var nhem_var shem_var;do fi done -# 1.c link in topo files -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.a ${DATA}/regional.grid.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.b ${DATA}/regional.grid.b -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.a ${DATA}/regional.depth.a -ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.b ${DATA}/regional.depth.b +# 1.c 1st try, don't link in MOM topo files +#ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.a ${DATA}/regional.grid.a +#ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.grid.b ${DATA}/regional.grid.b +#ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.a ${DATA}/regional.depth.a +#ln -f -s ${FIXrtofs}/${RUN}_${modID}.${inputgrid}.regional.depth.b ${DATA}/regional.depth.b # 2. Combine pre-qc and qc into one stream for ice and surface_obs/profile echo timecheck RTOFS_GLO_NCODA_QC start qc at $(date) diff --git a/sorc/hycom_tools_for_mom6.fd/.gitignore b/sorc/hycom_tools_for_mom6.fd/.gitignore new file mode 100644 index 00000000..e06b1478 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/.gitignore @@ -0,0 +1,5 @@ +*.o + +*.mod + +exec/* diff --git a/sorc/hycom_tools_for_mom6.fd/Makefile b/sorc/hycom_tools_for_mom6.fd/Makefile new file mode 100644 index 00000000..cebf9c1e --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/Makefile @@ -0,0 +1,70 @@ +# Makefile for compiling ncodaz_inc2mom6nc_glb_lyr + +# Add netCDF library paths and flags +NETCDF_PATH := $(NETCDF) +NETCDF_INC := $(NETCDF_PATH)/include +NETCDF_LIB := $(NETCDF_PATH)/lib + +# Compiler and flags +FC = ftn + +CRAY_INTEL_FFLAGS = -O1 -fpp -pc80 -traceback -WB -DCRAY +CRAY_INTEL_INTSIZE_4 = -i4 +CRAY_INTEL_REALSIZE_4 = -real-size 32 +CRAY_INTEL_SWAPIO = -convert big_endian -assume byterecl + +# Combined FFLAGS (Base flags for all files) +FFLAGS = $(CRAY_INTEL_FFLAGS) $(CRAY_INTEL_SWAPIO) \ + $(CRAY_INTEL_REALSIZE_4) $(CRAY_INTEL_INTSIZE_4) -I$(NETCDF_INC) + +# Linker flags - include netCDF libraries +LDFLAGS := -L$(NETCDF_LIB) -lnetcdf -lnetcdff + +# Source files (For reference) +SOURCES = ncodaz_inc2mom6nc_glb_lyr.f mod_mom6.F mod_za.F mod_ncio.f90 \ + mod_xc.F blkin.f zh.F getdat.f wtime.F + +# ALL objects must be listed here to satisfy the linker and pattern rules +#OBJECTS = ncodaz_inc2mom6nc_glb_lyr.o mod_mom6.o mod_za.o mod_ncio.o \ +# mod_xc.o blkin.o zh.o getdat.o wtime.o +OBJECTS = mod_xc.o mod_za.o mod_mom6.o mod_ncio.o blkin.o zh.o getdat.o wtime.o ncodaz_inc2mom6nc_glb_lyr.o + +EXECUTABLE = ncodaz_inc2mom6nc_glb_lyr.x + +# Default target +all: $(EXECUTABLE) + +# Link object files to create executable +$(EXECUTABLE): $(OBJECTS) + $(FC) $(FFLAGS) -o $@ $^ $(LDFLAGS) + +# --- Pattern Rules for different extensions --- + +# 1. Rule for legacy Fixed Form files (e.g., .f) +# This handles getdat.f, blkin.f, and ncodaz_inc2mom6nc_glb_lyr.f +%.o: %.f + $(FC) $(FFLAGS) -c $< + +# 2. Rule for Preprocessed / Free Form files (e.g., .F) +# Forces -free flag for mod_xc.F, mod_za.F, wtime.F, mod_mom6.F +%.o: %.F + $(FC) $(FFLAGS) -c $< +# $(FC) $(FFLAGS) -free -c $< + +# 3. Rule for modern f90 files (e.g., .f90 / .F90) +%.o: %.f90 + $(FC) $(FFLAGS) -free -c $< + +%.o: %.F90 + $(FC) $(FFLAGS) -free -c $< + +# Clean up object files and executable +clean: + rm -f $(OBJECTS) *.mod $(EXECUTABLE) + +# Remove only object files +clean-obj: + rm -f $(OBJECTS) + +# Phony targets +.PHONY: all clean clean-obj diff --git a/sorc/hycom_tools_for_mom6.fd/README.md b/sorc/hycom_tools_for_mom6.fd/README.md new file mode 100644 index 00000000..890709f9 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/README.md @@ -0,0 +1,2 @@ +# hycom_tools_for_RTOFS_v3 +Applications from HYCOM_tools that will be used in RTOFS version 3 with MOM6. diff --git a/sorc/hycom_tools_for_mom6.fd/bigrid.f b/sorc/hycom_tools_for_mom6.fd/bigrid.f new file mode 100644 index 00000000..d3bb99e7 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/bigrid.f @@ -0,0 +1,69 @@ + subroutine bigrid + use mod_mom6 ! HYCOM mom6 array interface + use mod_za ! HYCOM array I/O interface +c +c --- set mask arrays for irregular basin in c-grid configuration +c --- q,u,v,p are vorticity, u-velocity, v-velocity, and mass points, resp. +c --- 'depth' = basin depth array, zero values indicate land +c + write (lp,'(a,2i6)') + & 'bigrid called with ii,jj =', + & ii,jj + call zhflsh(lp) +c +c --- mass points are defined where water depth is greater than zero + do i=0,ii + do j=0,jj + if (depths(i,j).gt.0.) then + ip0(i,j)=1 + else + ip0(i,j)=0 + endif + enddo + enddo +c + do i=1,ii + do j=1,jj + ip(i,j)=ip0(i,j) + enddo + enddo +c +c --- u,v points are located halfway between any 2 adjoining mass points + do j=1,jj + do i=1,ii + if (ip0(i-1,j).gt.0.and.ip0(i,j).gt.0) then + iu(i,j)=1 + else + iu(i,j)=0 + endif + if (ip0(i,j-1).gt.0.and.ip0(i,j).gt.0) then + iv(i,j)=1 + else + iv(i,j)=0 + endif + enddo + enddo +c + do j=1,jj + do i=1,ii + iq(i,j)=0 +c +c --- 'interior' q points require water on all 4 sides. + if (min0(ip(i,j), ip(i-1,j), + & ip(i,j-1),ip(i-1,j-1)).gt.0) iq(i,j)=1 +c +c --- 'promontory' q points require water on 3 +c --- (or at least 2 diametrically opposed) sides + if ((ip(i ,j).gt.0.and.ip(i-1,j-1).gt.0).or. + & (ip(i-1,j).gt.0.and.ip(i ,j-1).gt.0) ) iq(i,j)=1 + enddo + enddo +c +* write(lp,'(a (/5x,70i1))') 'ip = ',(ip(i,1),i=1,ii) +* write(lp,'(a (/5x,70i1))') 'iu = ',(iu(i,1),i=1,ii) +* write(lp,'(a (/5x,70i1))') 'iv = ',(iv(i,1),i=1,ii) +* write(lp,'(a (/5x,70i1))') 'iq = ',(iq(i,1),i=1,ii) +* call zhflsh(lp) +c + return + end diff --git a/sorc/hycom_tools_for_mom6.fd/blkin.f b/sorc/hycom_tools_for_mom6.fd/blkin.f new file mode 100644 index 00000000..9bbfb2eb --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/blkin.f @@ -0,0 +1,442 @@ + subroutine blkind(dvar,cvar,cfmt) + implicit none +c + double precision dvar + character cvar*6,cfmt*(*) +c + integer lp + common/linepr/lp +c +c read in one d.p. value from stdin +c + character*6 cvarin +c + read(*,*) dvar,cvarin + write(lp,cfmt) cvarin,dvar + call flush(lp) +c + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkind - input ',cvarin, + + ' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + end + subroutine blkind2(dvar,nvar,cvar1,cfmt1,cvar2,cfmt2) + implicit none +c + double precision dvar + integer nvar + character cvar1*6,cvar2*6,cfmt1*(*),cfmt2*(*) +c + integer lp + common/linepr/lp +c +c read in one d.p. value from stdin, +c identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) +c + character*6 cvarin +c + read(*,*) dvar,cvarin +c + if (cvar1.eq.cvarin) then + nvar = 1 + write(lp,cfmt1) cvarin,dvar + call flush(lp) + elseif (cvar2.eq.cvarin) then + nvar = 2 + write(lp,cfmt2) cvarin,dvar + call flush(lp) + else + write(lp,*) + write(lp,*) 'error in blkind2 - input ',cvarin, + + ' but should be ',cvar1,' or ',cvar2 + write(lp,*) + call flush(lp) + stop + endif + return + end + subroutine blkinr(rvar,cvar,cfmt) + implicit none +c + real rvar + character cvar*6,cfmt*(*) +c + integer lp + common/linepr/lp +c +c read in one real value from stdin +c + character*6 cvarin +c + read(*,*) rvar,cvarin + write(lp,cfmt) cvarin,rvar + call flush(lp) +c + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkinr - input ',cvarin, + + ' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + end + subroutine blkinr2(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2) + implicit none +c + real rvar + integer nvar + character cvar1*6,cvar2*6,cfmt1*(*),cfmt2*(*) +c + integer lp + common/linepr/lp +c +c read in one real value from stdin, +c identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) +c + character*6 cvarin +c + read(*,*) rvar,cvarin +c + if (cvar1.eq.cvarin) then + nvar = 1 + write(lp,cfmt1) cvarin,rvar + call flush(lp) + elseif (cvar2.eq.cvarin) then + nvar = 2 + write(lp,cfmt2) cvarin,rvar + call flush(lp) + else + write(lp,*) + write(lp,*) 'error in blkinr2 - input ',cvarin, + + ' but should be ',cvar1,' or ',cvar2 + write(lp,*) + call flush(lp) + stop + endif + return + end + subroutine blkinr3(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3) + implicit none +c + real rvar + integer nvar + character cvar1*6,cvar2*6,cvar3*6,cfmt1*(*),cfmt2*(*),cfmt3*(*) +c + integer lp + common/linepr/lp +c +c read in one of three possible real values from stdin, +c identified as either cvar1 (return nvar=1) or +c cvar2 (return nvar=2) or +c cvar3 (return nvar=3) +c + call blkinr9(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3, + & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', + & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', + & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', + & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', + & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', + & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")') + return + end + subroutine blkinr9(rvar,nvar, + & cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3, + & cvar4,cfmt4,cvar5,cfmt5,cvar6,cfmt6, + & cvar7,cfmt7,cvar8,cfmt8,cvar9,cfmt9) + implicit none +c + real rvar + integer nvar + character cvar1*6,cvar2*6,cvar3*6,cfmt1*(*),cfmt2*(*),cfmt3*(*) + character cvar4*6,cvar5*6,cvar6*6,cfmt4*(*),cfmt5*(*),cfmt6*(*) + character cvar7*6,cvar8*6,cvar9*6,cfmt7*(*),cfmt8*(*),cfmt9*(*) +c + integer lp + common/linepr/lp +c +c read in one of nine possible real values from stdin, +c identified as either cvar1 (return nvar=1) or +c cvar2 (return nvar=2) or +c ... +c cvar9 (return nvar=9) +c + character*6 cvarin +c + read(*,*) rvar,cvarin +c + if (cvar1.eq.cvarin) then + nvar = 1 + write(lp,cfmt1) cvarin,rvar + call flush(lp) + elseif (cvar2.eq.cvarin) then + nvar = 2 + write(lp,cfmt2) cvarin,rvar + call flush(lp) + elseif (cvar3.eq.cvarin) then + nvar = 3 + write(lp,cfmt3) cvarin,rvar + call flush(lp) + elseif (cvar4.eq.cvarin) then + nvar = 4 + write(lp,cfmt4) cvarin,rvar + call flush(lp) + elseif (cvar5.eq.cvarin) then + nvar = 5 + write(lp,cfmt5) cvarin,rvar + call flush(lp) + elseif (cvar6.eq.cvarin) then + nvar = 6 + write(lp,cfmt6) cvarin,rvar + call flush(lp) + elseif (cvar7.eq.cvarin) then + nvar = 7 + write(lp,cfmt7) cvarin,rvar + call flush(lp) + elseif (cvar8.eq.cvarin) then + nvar = 8 + write(lp,cfmt8) cvarin,rvar + call flush(lp) + elseif (cvar9.eq.cvarin) then + nvar = 9 + write(lp,cfmt9) cvarin,rvar + call flush(lp) + else + write(lp,*) + write(lp,*) 'error in blkinr9 - input ',cvarin, + + ' but should be:' + write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or' + write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or' + write(lp,*) cvar7,' or ',cvar8,' or ',cvar9 + write(lp,*) + call flush(lp) + stop + endif + return + end + subroutine blkini(ivar,cvar) + implicit none +c + integer ivar + character*6 cvar +c + integer lp + common/linepr/lp +c +c read in one integer value from stdin +c + character*6 cvarin +c + read(*,*) ivar,cvarin + write(lp,6000) cvarin,ivar + call flush(lp) +c + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkini - input ',cvarin, + + ' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkini: ',a6,' =',i6) + end + subroutine blkini2(ivar,nvar,cvar1,cvar2) + implicit none +c + integer ivar,nvar + character*6 cvar1,cvar2 +c + integer lp + common/linepr/lp +c +c read in one integer value from stdin +c identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) +c + character*6 cvarin +c + read(*,*) ivar,cvarin + write(lp,6000) cvarin,ivar + call flush(lp) +c + if (cvarin.eq.cvar1) then + nvar = 1 + elseif (cvarin.eq.cvar2) then + nvar = 2 + else + write(lp,*) + write(lp,*) 'error in blkini2 - input ',cvarin, + + ' but should be ',cvar1,' or ',cvar2 + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkini: ',a6,' =',i6) + end + subroutine blkini3(ivar,nvar,cvar1,cvar2,cvar3) + implicit none +c + integer ivar,nvar + character*6 cvar1,cvar2,cvar3 +c + integer lp + common/linepr/lp +c +c read in one of three possible integer values from stdin +c identified as either cvar1 (return nvar=1) or +c cvar2 (return nvar=2) or +c cvar3 (return nvar=3) +c + call blkini9(ivar,nvar,cvar1,cvar2,cvar3, + & 'XXXXXX','XXXXXX','XXXXXX', + & 'XXXXXX','XXXXXX','XXXXXX') + return + end + subroutine blkini4(ivar,nvar,cvar1,cvar2,cvar3,cvar4) + implicit none +c + integer ivar,nvar + character*6 cvar1,cvar2,cvar3,cvar4 +c + integer lp + common/linepr/lp +c +c read in one of four possible integer values from stdin +c identified as either cvar1 (return nvar=1) or +c cvar2 (return nvar=2) or +c cvar3 (return nvar=3) or +c cvar4 (return nvar=4) +c + call blkini9(ivar,nvar,cvar1, cvar2, cvar3, + & cvar4, 'XXXXXX','XXXXXX', + & 'XXXXXX','XXXXXX','XXXXXX') + return + end + subroutine blkini9(ivar,nvar,cvar1,cvar2,cvar3, + & cvar4,cvar5,cvar6, + & cvar7,cvar8,cvar9) + implicit none +c + integer ivar,nvar + character*6 cvar1,cvar2,cvar3 + character*6 cvar4,cvar5,cvar6 + character*6 cvar7,cvar8,cvar9 +c + integer lp + common/linepr/lp +c +c read in one of nine possible integer values from stdin +c identified as either cvar1 (return nvar=1) or +c cvar2 (return nvar=2) or +c ... +c cvar9 (return nvar=9) +c cvar3 (return nvar=3) or +c + character*6 cvarin +c + read(*,*) ivar,cvarin + write(lp,6000) cvarin,ivar + call flush(lp) +c + if (cvarin.eq.cvar1) then + nvar = 1 + elseif (cvarin.eq.cvar2) then + nvar = 2 + elseif (cvarin.eq.cvar3) then + nvar = 3 + elseif (cvar4.eq.cvarin) then + nvar = 4 + elseif (cvar5.eq.cvarin) then + nvar = 5 + elseif (cvar6.eq.cvarin) then + nvar = 6 + elseif (cvar7.eq.cvarin) then + nvar = 7 + elseif (cvar8.eq.cvarin) then + nvar = 8 + elseif (cvar9.eq.cvarin) then + nvar = 9 + else + write(lp,*) + write(lp,*) 'error in blkini9 - input ',cvarin, + + ' but should be:' + write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or' + write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or' + write(lp,*) cvar7,' or ',cvar8,' or ',cvar9 + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkini: ',a6,' =',i6) + end + subroutine blkinl(lvar,cvar) + implicit none +c + logical lvar + character*6 cvar +c + integer lp + common/linepr/lp +c +c read in one logical value from stdin +c due to a SGI bug for logical I/O: read in an integer 0=F,1=T +c + character*6 cvarin + integer ivar +c + read(*,*) ivar,cvarin + lvar = ivar .ne. 0 + write(lp,6000) cvarin,lvar + call flush(lp) +c + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkinr - input ',cvarin, + + ' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkinl: ',a6,' =',l6) + end + subroutine blkinl_99(lvar,cvar) + implicit none +c + logical lvar + character*6 cvar +c + integer lp + common/linepr/lp +c +c read in one logical value from unit 99 +c due to a SGI bug for logical I/O: read in an integer 0=F,1=T +c + character*6 cvarin + integer ivar +c + read(99,*) ivar,cvarin + lvar = ivar .ne. 0 + write(lp,6000) cvarin,lvar + call flush(lp) +c + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkinr - input ',cvarin, + + ' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkinl: ',a6,' =',l6) + end diff --git a/sorc/hycom_tools_for_mom6.fd/build_all.sh b/sorc/hycom_tools_for_mom6.fd/build_all.sh new file mode 100755 index 00000000..9dbef1ef --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/build_all.sh @@ -0,0 +1,28 @@ +#!/bin/sh + +BASE=`pwd` +dir_mods="$(dirname ${BASE})" +dir_mod0="$(dirname ${dir_mods})" + +echo " " +echo "Load modules listed at: " +echo ${dir_mod0}"/versions/build.ver" +source ${BASE}/load_modules.sh ${dir_mod0} + +make all + +# move executables +mkdir -p ${BASE}/exec +mv *.x ${BASE}/exec + +#clean: $(SUBDIRS) +# for dir in $(SUBDIRS); do \ +# ( cd $$dir; echo "Making $@ in `pwd`" ; \ +# make $@) ; \ +# done + +#install: $(SUBDIRS) +# for dir in $(SUBDIRS); do \ +# (cd $$dir; cp -p bin/* ../../../exec); \ +# done + diff --git a/sorc/hycom_tools_for_mom6.fd/dum_gks.f b/sorc/hycom_tools_for_mom6.fd/dum_gks.f new file mode 100644 index 00000000..7438c50f --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/dum_gks.f @@ -0,0 +1,6 @@ + subroutine clsgks +c +c --- dummy version for no graphics cases. +c + return + end diff --git a/sorc/hycom_tools_for_mom6.fd/forday.f b/sorc/hycom_tools_for_mom6.fd/forday.f new file mode 100644 index 00000000..421f5c3f --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/forday.f @@ -0,0 +1,147 @@ + subroutine fordate(dtime,yrflag, iyear,month,iday,ihour) + implicit none +c + double precision dtime + integer yrflag, iyear,month,iday,ihour +c +c --- converts model day to "calendar" date (year,month,day,hour). +c + integer jday,k,m +c + integer month0(13,3) + data month0 / 1, 31, 61, 91, 121, 151, 181, + + 211, 241, 271, 301, 331, 361, + + 1, 32, 60, 91, 121, 152, 182, + + 213, 244, 274, 305, 335, 366, + + 1, 32, 61, 92, 122, 153, 183, + + 214, 245, 275, 306, 336, 367 / +c + call forday(dtime,yrflag, iyear,jday,ihour) +c + if (yrflag.eq.3) then + if (mod(iyear,4).eq.0) then + k = 3 + else + k = 2 + endif + elseif (yrflag.eq.0) then + k = 1 + else + k = 3 + endif + do m= 1,12 + if (jday.ge.month0(m, k) .and. + + jday.lt.month0(m+1,k) ) then + month = m + iday = jday - month0(m,k) + 1 + endif + enddo + return + end + subroutine forday(dtime,yrflag, iyear,iday,ihour) + implicit none +c + double precision dtime + integer yrflag, iyear,iday,ihour +c +c --- converts model day to "calendar" date (year,julian-day,hour). +c + integer lp + common/linepr/ lp + save /linepr/ +c + double precision dtim1,day + integer iyr,nleap +c + if (yrflag.eq.0) then +c --- 360 days per model year, starting Jan 16 + iyear = int((dtime+15.001d0)/360.d0) + 1 + iday = mod( dtime+15.001d0 ,360.d0) + 1 + ihour = (mod( dtime+15.001d0 ,360.d0) + 1.d0 - iday)*24.d0 +c + elseif (yrflag.eq.1) then +c --- 366 days per model year, starting Jan 16 + iyear = int((dtime+15.001d0)/366.d0) + 1 + iday = mod( dtime+15.001d0 ,366.d0) + 1 + ihour = (mod( dtime+15.001d0 ,366.d0) + 1.d0 - iday)*24.d0 +c + elseif (yrflag.eq.2) then +c --- 366 days per model year, starting Jan 01 + iyear = int((dtime+ 0.001d0)/366.d0) + 1 + iday = mod( dtime+ 0.001d0 ,366.d0) + 1 + ihour = (mod( dtime+ 0.001d0 ,366.d0) + 1.d0 - iday)*24.d0 +c + elseif (yrflag.eq.3) then +c --- model day is calendar days since 01/01/1901 + iyr = (dtime-1.d0)/365.25d0 + nleap = iyr/4 + dtim1 = 365.d0*iyr + nleap + 1.d0 + day = dtime - dtim1 + 1.d0 + if (dtim1.gt.dtime) then + iyr = iyr - 1 + elseif (day.ge.367.d0) then + iyr = iyr + 1 + elseif (day.ge.366.d0 .and. mod(iyr,4).ne.3) then + iyr = iyr + 1 + endif + nleap = iyr/4 + dtim1 = 365.d0*iyr + nleap + 1.d0 +c + iyear = 1901 + iyr + iday = dtime - dtim1 + 1 + ihour = (dtime - dtim1 + 1.d0 - iday)*24.d0 +c + else + write(lp,*) + write(lp,*) 'error in forday - unsupported yrflag value' + write(lp,*) + call flush(lp) + stop '(forday)' + endif + return + end + + SUBROUTINE DATE2WNDAY(WDAY, IDATEC,ITIMEC) + IMPLICIT NONE + INTEGER IDATEC,ITIMEC + DOUBLE PRECISION WDAY +C +C********** +C* +C 1) CONVERT DATE INTO 'FLUX DAY'. +C +C 2) THE 'FLUX DAY' IS THE NUMBER OF DAYS SINCE 001/1901 (WHICH IS +C FLUX DAY 1.0). +C FOR EXAMPLE: +C A) IYR=1901,MON=1,IDY=1, REPRESENTS 0000Z HRS ON 01/01/1901 +C SO WDAY WOULD BE 1.0. +C A) IYR=1901,MON=1,IDY=2, REPRESENTS 0000Z HRS ON 02/01/1901 +C SO WDAY WOULD BE 2.0. +C YEAR MUST BE NO LESS THAN 1901.0, AND NO GREATER THAN 2099.0. +C NOTE THAT YEAR 2000 IS A LEAP YEAR (BUT 1900 AND 2100 ARE NOT). +C +C 3) ALAN J. WALLCRAFT, NAVAL RESEARCH LABORATORY, JULY 2002. +C* +C********** +C + INTEGER NLEAP,IYR,MON,IDY,IHR +C + INTEGER MONTH(13) + DATA MONTH / 0, 31, 59, 90, 120, 151, 181, + + 212, 243, 273, 304, 334, 365 / +C + IYR = IDATEC/10000 + MON = (IDATEC-IYR*10000)/100 + IDY = MOD(IDATEC,100) + IHR = ITIMEC/1000000 +C +C FIND THE RIGHT YEAR. +C + NLEAP = (IYR-1901)/4 + WDAY = 365.0*(IYR-1901) + NLEAP + MONTH(MON) + IDY + IHR/24.0 + IF (MOD(IYR,4).EQ.0 .AND. MON.GT.2) THEN + WDAY = WDAY + 1.0 + ENDIF + RETURN +C END OF DATE2WNDAY + END diff --git a/sorc/hycom_tools_for_mom6.fd/getdat.f b/sorc/hycom_tools_for_mom6.fd/getdat.f new file mode 100644 index 00000000..b8505af2 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/getdat.f @@ -0,0 +1,1441 @@ + subroutine rd_out2nc(n,m,irec, + & t,time3, + & name_t,flnm_t) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + character*(*) name_t,flnm_t + integer n,m,irec + double precision time3(3) + real t(n,m) +c +c subroutine to read a MOM6 scalar 2-d field. +c +c subroutine arguments: +c n,m = horizontal grid dimensions. +c irec = time record to input (0-> input record "time") +c unchanged on output if >0, set to read record otherwise +c flnm_t = filename of the netCDF input. +c name_t = name of the netCDF field. +c +c t = field +c time3 = days since 1901-01-01 00:00:00 +c time3(1) = start time interval +c time3(2) = end time interval +c time3(3) = mid time interval +c output if irec>0, input if irec=0 +c + character*(NF90_MAX_NAME) :: cD + character*(256) :: flnm_p + logical :: lexist + integer :: ncFID,ncVID,ncNID,ncSTATUS + integer :: ncDIDs(nf90_max_var_dims) + integer :: i,k,nc,nc1st,nclst + integer :: np,mp,nx(4),ny(4),tto + real, allocatable :: p(:,:) + double precision, allocatable :: t_rec(:) + double precision :: tbnds(2) +c + write(6,*) "irec = ",irec +c + i = len_trim(flnm_t) + if (flnm_t(i-2:i).eq.'.nc') then +c --- single netCDF file + call nchek('nf90_open-T', + & nf90_open(trim(flnm_t), nf90_nowrite, ncFID)) + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(2) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + write(6,*) "name = ",trim(name_t) + write(6,*) "time3 = ",time3 + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, t(:,:), + & (/ 1,1,irec /) )) + call nchek("nf90_close", + & nf90_close(ncFID)) + else +c --- multiple netCDF files + read(flnm_t(i-8:i-5),*) nc1st + read(flnm_t(i-3:i) ,*) nclst + do nc= nc1st, nclst +* write (6,*) 'nc = ',nc + flnm_p = flnm_t(1:i-9) + write(flnm_p(i-8:i-5),'(i4.4)') nc + lexist = nf90_open(trim(flnm_p), nf90_nowrite, ncFID) + & .eq. nf90_noerr + if (lexist) then + write (6,'(2a)') ' input MOM6 file: ',trim(flnm_p) + else + write (6,'(2a)') ' skip MOM6 file: ',trim(flnm_p) + cycle + endif + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), + & len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(2) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + if (nc.eq.nc1st) then + write(6,*) "name = ",trim(name_t) + write(6,*) "time3 = ",time3 + endif +c + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", nx(:))) +* write(6,*) 'nx = ',nx +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", ny(:))) +* write(6,*) 'ny = ',ny +c +c --- symetric arrays may have an extra row or column +c + nx(4) = min(nx(4),n) + ny(4) = min(ny(4),m) +c +c --- symetric arrays may have an extra row or column +c + nx(4) = min(nx(4),n) + ny(4) = min(ny(4),m) +c + np = nx(4)-nx(3)+1 + mp = ny(4)-ny(3)+1 +* write(6,*) 'np,mp = ',np,mp + allocate( p(np,mp) ) + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, p(:,:), + & (/ 1,1,irec /) )) + t(nx(3):nx(4),ny(3):ny(4)) = p(:,:) + deallocate( p ) + call nchek("nf90_close", + & nf90_close(ncFID)) + enddo !nc + endif +c + return + end + + subroutine rd_out2nc8(n,m,irec, + & t,time3, + & name_t,flnm_t) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + character*(*) name_t,flnm_t + integer n,m,irec + double precision time3(3) + double precision t(n,m) +c +c subroutine to read a MOM6 scalar double 2-d field. +c +c subroutine arguments: +c n,m = horizontal grid dimensions. +c irec = time record to input (0-> input record "time") +c unchanged on output if >0, set to read record otherwise +c flnm_t = filename of the netCDF input. +c name_t = name of the netCDF field. +c +c t = field (real*8) +c time3 = days since 1901-01-01 00:00:00 +c time3(1) = start time interval +c time3(2) = end time interval +c time3(3) = mid time interval +c output if irec>0, input if irec=0 +c + character*(NF90_MAX_NAME) :: cD + character*(256) :: flnm_p + logical :: lexist + integer :: ncFID,ncVID,ncNID,ncSTATUS + integer :: ncDIDs(nf90_max_var_dims) + integer :: i,k,nc,nc1st,nclst + integer :: np,mp,nx(4),ny(4),tto + real, allocatable :: p(:,:) + double precision, allocatable :: t_rec(:) + double precision :: tbnds(2) +c + write(6,*) "irec = ",irec +c + i = len_trim(flnm_t) + if (flnm_t(i-2:i).eq.'.nc') then +c --- single netCDF file + call nchek('nf90_open-T', + & nf90_open(trim(flnm_t), nf90_nowrite, ncFID)) + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(2) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + write(6,*) "name = ",trim(name_t) + write(6,*) "time3 = ",time3 + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, t(:,:), + & (/ 1,1,irec /) )) + call nchek("nf90_close", + & nf90_close(ncFID)) + else +c --- multiple netCDF files + read(flnm_t(i-8:i-5),*) nc1st + read(flnm_t(i-3:i) ,*) nclst + do nc= nc1st, nclst +* write (6,*) 'nc = ',nc + flnm_p = flnm_t(1:i-9) + write(flnm_p(i-8:i-5),'(i4.4)') nc + lexist = nf90_open(trim(flnm_p), nf90_nowrite, ncFID) + & .eq. nf90_noerr + if (lexist) then + write (6,'(2a)') ' input MOM6 file: ',trim(flnm_p) + else + write (6,'(2a)') ' skip MOM6 file: ',trim(flnm_p) + cycle + endif + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), + & len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(2) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + if (nc.eq.nc1st) then + write(6,*) "name = ",trim(name_t) + write(6,*) "time3 = ",time3 + endif +c + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", nx(:))) +* write(6,*) 'nx = ',nx +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", ny(:))) +* write(6,*) 'ny = ',ny +c +c --- symetric arrays may have an extra row or column +c + nx(4) = min(nx(4),n) + ny(4) = min(ny(4),m) +c + np = nx(4)-nx(3)+1 + mp = ny(4)-ny(3)+1 +* write(6,*) 'np,mp = ',np,mp + allocate( p(np,mp) ) + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, p(:,:), + & (/ 1,1,irec /) )) + t(nx(3):nx(4),ny(3):ny(4)) = p(:,:) + deallocate( p ) + call nchek("nf90_close", + & nf90_close(ncFID)) + enddo !nc + endif +c + return + end + + subroutine rd_out3nc(n,m,l,irec, + & t,time3, + & name_t,flnm_t) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + character*(*) name_t,flnm_t + integer n,m,l,irec + double precision time3(3) + real t(n,m,l) +c +c subroutine to read a MOM6 3-d field. +c +c subroutine arguments: +c n,m = horizontal grid dimensions. +c l = vertical grid dimension. +c irec = time record to input (0-> input record "time") +c unchanged on output if >0, set to read record otherwise +c flnm_t = filename of the netCDF input. +c name_t = name of the netCDF field. +c +c t = field +c time3 = days since 1901-01-01 00:00:00 +c time3(1) = start time interval +c time3(2) = end time interval +c time3(3) = mid time interval +c output if irec>0, input if irec=0 +c + character*(NF90_MAX_NAME) :: cD + character*(256) :: flnm_p + logical :: lexist + integer :: ncFID,ncVID,ncNID,ncSTATUS + integer :: ncDIDs(nf90_max_var_dims) + integer :: i,k,nc,nc1st,nclst + integer :: np,mp,nx(4),ny(4),tto + real, allocatable :: p(:,:,:) + double precision, allocatable :: t_rec(:) + double precision :: tbnds(2) +c + write(6,*) "irec = ",irec +c + i = len_trim(flnm_t) + if (flnm_t(i-2:i).eq.'.nc') then +c --- single netCDF file + call nchek('nf90_open-T', + & nf90_open(trim(flnm_t), nf90_nowrite, ncFID)) + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(1) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + write(6,*) "name = ",trim(name_t) + write(6,*) "time = ",time3(3) + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, t(:,:,:), + & (/ 1,1,1,irec /) )) + call nchek("nf90_close", + & nf90_close(ncFID)) + else +c --- multiple netCDF files + read(flnm_t(i-8:i-5),*) nc1st + read(flnm_t(i-3:i) ,*) nclst + do nc= nc1st, nclst +* write (6,*) 'nc = ',nc + flnm_p = flnm_t(1:i-9) + write(flnm_p(i-8:i-5),'(i4.4)') nc + lexist = nf90_open(trim(flnm_p), nf90_nowrite, ncFID) + & .eq. nf90_noerr + if (lexist) then + write (6,'(2a)') ' input MOM6 file: ',trim(flnm_p) + else + write (6,'(2a)') ' skip MOM6 file: ',trim(flnm_p) + cycle + endif + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), + & len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(1) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + if (nc.eq.nc1st) then + write(6,*) "name = ",trim(name_t) + write(6,*) "time3 = ",time3 + endif +c + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", nx(:))) +* write(6,*) 'nx = ',nx +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", ny(:))) +* write(6,*) 'ny = ',ny +c +c --- symetric arrays may have an extra row or column +c + nx(4) = min(nx(4),n) + ny(4) = min(ny(4),m) +c + np = nx(4)-nx(3)+1 + mp = ny(4)-ny(3)+1 +* write(6,*) 'np,mp = ',np,mp + allocate( p(np,mp,l) ) + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, + & p(:,:,:), + & (/ 1,1,1,irec /) )) + t(nx(3):nx(4),ny(3):ny(4),1:l) = p(:,:,:) + deallocate( p ) + call nchek("nf90_close", + & nf90_close(ncFID)) + enddo !nc + endif +c + return + end + + subroutine rd_out3nc8(n,m,l,irec, + & t,time3, + & name_t,flnm_t) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + character*(*) name_t,flnm_t + integer n,m,l,irec + double precision time3(3) + double precision t(n,m,l) +c +c subroutine to read a MOM6 3-d double field. +c +c subroutine arguments: +c n,m = horizontal grid dimensions. +c l = vertical grid dimension. +c irec = time record to input (0-> input record "time") +c unchanged on output if >0, set to read record otherwise +c flnm_t = filename of the netCDF input. +c name_t = name of the netCDF field. +c +c t = field (real*8) +c time3 = days since 1901-01-01 00:00:00 +c time3(1) = start time interval +c time3(2) = end time interval +c time3(3) = mid time interval +c output if irec>0, input if irec=0 +c + character*(NF90_MAX_NAME) :: cD + character*(256) :: flnm_p + logical :: lexist + integer :: ncFID,ncVID,ncNID,ncSTATUS + integer :: ncDIDs(nf90_max_var_dims) + integer :: i,k,nc,nc1st,nclst + integer :: np,mp,nx(4),ny(4),tto + real, allocatable :: p(:,:,:) + double precision, allocatable :: t_rec(:) + double precision :: tbnds(2) +c + write(6,*) "irec = ",irec +c + i = len_trim(flnm_t) + if (flnm_t(i-2:i).eq.'.nc') then +c --- single netCDF file + call nchek('nf90_open-T', + & nf90_open(trim(flnm_t), nf90_nowrite, ncFID)) + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(1) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + write(6,*) "name = ",trim(name_t) + write(6,*) "time = ",time3(3) + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, t(:,:,:), + & (/ 1,1,1,irec /) )) + call nchek("nf90_close", + & nf90_close(ncFID)) + else +c --- multiple netCDF files + read(flnm_t(i-8:i-5),*) nc1st + read(flnm_t(i-3:i) ,*) nclst + do nc= nc1st, nclst +* write (6,*) 'nc = ',nc + flnm_p = flnm_t(1:i-9) + write(flnm_p(i-8:i-5),'(i4.4)') nc + lexist = nf90_open(trim(flnm_p), nf90_nowrite, ncFID) + & .eq. nf90_noerr + if (lexist) then + write (6,'(2a)') ' input MOM6 file: ',trim(flnm_p) + else + write (6,'(2a)') ' skip MOM6 file: ',trim(flnm_p) + cycle + endif + call nchek('nf90_inq_varid-Time', + & nf90_inq_varid(ncFID,'Time', ncVID)) + if (irec.eq.0) then + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(1:1))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), + & len=tto)) + allocate( t_rec(tto) ) + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, t_rec(:), + & (/ 1 /) )) + do k= 1,tto + if (abs(t_rec(k)-time3(3)).lt.0.02) then + irec = k + write(6,*) "irec = ",irec + exit + endif + enddo !k + if (irec.eq.0) then + write(6,*) + write(6,*) 'error - requested time not in time dimension' + write(6,*) 'time = ',time3(3) + write(6,*) 't_rec = ',t_rec(:) + write(6,*) + stop + endif !irec==0 + deallocate( t_rec ) + else + call nchek('nf90_get_var-Time', + & nf90_get_var( ncFID, ncVID, time3(3), + & (/ irec /) )) + endif + ncSTATUS = nf90_inq_varid(ncFID,'Time_bnds', ncVID) + if (ncSTATUS.eq.NF90_NOERR) then +c --- time-mean fields + call nchek('nf90_get_var-Time_bnds', + & nf90_get_var( ncFID, ncVID, tbnds(1:2), + & (/ 1,irec /) )) + time3(1) = tbnds(1) + time3(2) = tbnds(1) + else +c --- instantaneous fields + time3(1) = time3(3) + time3(2) = time3(3) + endif + if (nc.eq.nc1st) then + write(6,*) "name = ",trim(name_t) + write(6,*) "time3 = ",time3 + endif +c + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", nx(:))) +* write(6,*) 'nx = ',nx +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", ny(:))) +* write(6,*) 'ny = ',ny +c +c --- symetric arrays may have an extra row or column +c + nx(4) = min(nx(4),n) + ny(4) = min(ny(4),m) +c + np = nx(4)-nx(3)+1 + mp = ny(4)-ny(3)+1 +* write(6,*) 'np,mp = ',np,mp + allocate( p(np,mp,l) ) + call nchek('nf90_inq_varid-'//trim(name_t), + & nf90_inq_varid(ncFID,trim(name_t), ncVID)) + call nchek('nf90_get_var-'//trim(name_t), + & nf90_get_var( ncFID, ncVID, + & p(:,:,:), + & (/ 1,1,1,irec /) )) + t(nx(3):nx(4),ny(3):ny(4),1:l) = p(:,:,:) + deallocate( p ) + call nchek("nf90_close", + & nf90_close(ncFID)) + enddo !nc + endif +c + return + end + + subroutine rd_dimen(xto,yto,zto,tto, cfile,cname) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + integer xto,yto,zto,tto + character*(*) cfile,cname +c +c subroutine to read model dimensions +c xto,yto= horizontal dimensions of entire grid +c zto = total number of vertical layers +c tto = total number of time samples +c + character*(NF90_MAX_NAME) :: cD + character*(256) :: flnm_p + logical :: lexist + integer i,ncFID,ncDID,ncVID,ncNID,ncDIDs(nf90_max_var_dims) + integer nc,nc1st,nclst + integer nx(4),ny(4) +c + i = len_trim(cfile) + if (cfile(i-2:i).eq.'.nc') then + call nchek('nf90_open-rd_dimen', + & nf90_open(trim(cfile), nf90_nowrite, ncFID)) +c + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cname), ncVID)) +c + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) +c + if (ncNID.ne.4) then + write(6,'(/ 3a /)') + & 'error - variable ',trim(cname),' does not have 4 dimensions' + stop + endif +c + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) +c + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), len=xto)) +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), len=yto)) +c + call nchek('nf90_inquire_dimension-3', + & nf90_inquire_dimension(ncFID, ncDIDs(3), len=zto)) +c + call nchek('nf90_inquire_dimension-4', + & nf90_inquire_dimension(ncFID, ncDIDs(4), len=tto)) +c + call nchek("nf90_close", + & nf90_close(ncFID)) + else +c --- multiple netCDF files + read(cfile(i-8:i-5),*) nc1st + read(cfile(i-3:i) ,*) nclst + do nc= nc1st, nclst + write (6,*) 'nc = ',nc + flnm_p = cfile(1:i-9) + write(flnm_p(i-8:i-5),'(i4.4)') nc + lexist = nf90_open(trim(flnm_p), nf90_nowrite, ncFID) + & .eq. nf90_noerr + if (.not.lexist) then + if (nc.eq.nclst) then + write(6,'(/ 3a /)') + & 'error - no input files (',trim(cfile),')' + stop + endif + cycle !nc loop + endif + write(6,'(3a)') "nf90_open(",trim(flnm_p),"," + call nchek('nf90_open-rd_dimen-XXXX', + & nf90_open(trim(flnm_p), nf90_nowrite, ncFID)) +c + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cname), ncVID)) +c + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) +c + if (ncNID.ne.4) then + write(6,'(/ 3a /)') + & 'error - variable ',trim(cname),' does not have 4 dimensions' + stop + endif +c + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) +c + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", nx(:))) + write(6,*) 'nx = ',nx + xto = nx(2) +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", ny(:))) + write(6,*) 'ny = ',ny + yto = ny(2) +c + call nchek('nf90_inquire_dimension-3', + & nf90_inquire_dimension(ncFID, ncDIDs(3), len=zto)) +c + call nchek('nf90_inquire_dimension-4', + & nf90_inquire_dimension(ncFID, ncDIDs(4), len=tto)) +c + call nchek("nf90_close", + & nf90_close(ncFID)) +c --- found an input file + exit !nc loop + enddo !nc + endif + return + end + + subroutine rd_dimen2(xto,yto,tto, cfile,cname) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + integer xto,yto,tto + character*(*) cfile,cname +c +c subroutine to read model 2-D dimensions +c xto,yto= horizontal dimensions of entire grid +c tto = total number of time samples +c + character*(NF90_MAX_NAME) :: cD + integer i,ncFID,ncDID,ncVID,ncNID,ncDIDs(nf90_max_var_dims) + integer nx(4),ny(4) +c + i = len_trim(cfile) + if (cfile(i-2:i).eq.'.nc') then + call nchek('nf90_open-rd_dimen2', + & nf90_open(trim(cfile), nf90_nowrite, ncFID)) +c + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cname), ncVID)) +c + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) +c + if (ncNID.ne.3) then + write(6,'(/ 3a /)') + & 'error - variable ',trim(cname),' does not have 3 dimensions' + stop + endif +c + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) +c + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), len=xto)) +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), len=yto)) +c + call nchek('nf90_inquire_dimension-3', + & nf90_inquire_dimension(ncFID, ncDIDs(3), len=tto)) +c + call nchek("nf90_close", + & nf90_close(ncFID)) + else + i = len_trim(cfile) +* write(6,'(3a)') "nf90_open(",cfile(1:i-5),"," + call nchek('nf90_open-rd_dimen2-0000', + & nf90_open(cfile(1:i-5), nf90_nowrite, ncFID)) +c + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cname), ncVID)) +c + call nchek('nf90_inq_variable(ndims)', + & nf90_inquire_variable(ncFID, ncVID, ndims=ncNID)) +c + if (ncNID.ne.3) then + write(6,'(/ 3a /)') + & 'error - variable ',trim(cname),' does not have 3 dimensions' + stop + endif +c + call nchek('nf90_inq_variable(dimids)', + & nf90_inquire_variable(ncFID, ncVID, + & dimids=ncDIDs(:ncNID))) +c + call nchek('nf90_inquire_dimension-1', + & nf90_inquire_dimension(ncFID, ncDIDs(1), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", nx(:))) +* write(6,*) 'nx = ',nx + xto = nx(2) +c + call nchek('nf90_inquire_dimension-2', + & nf90_inquire_dimension(ncFID, ncDIDs(2), name=cD)) + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cD), ncVID)) + call nchek('nf90_get_att', + & nf90_get_att(ncFID, ncVID, + & "domain_decomposition", ny(:))) +* write(6,*) 'ny = ',ny + yto = ny(2) +c + call nchek('nf90_inquire_dimension-3', + & nf90_inquire_dimension(ncFID, ncDIDs(3), len=tto)) +c + call nchek("nf90_close", + & nf90_close(ncFID)) + endif + return + end + + subroutine rd_missing(misval, cfile,cname) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + real misval + character*(*) cfile,cname +c +c subroutine to read MOM6 missing_value +c misval = missing_value +c + character*(NF90_MAX_NAME) :: cD + integer i,ncFID,ncDID,ncVID,ncNID,ncDIDs(nf90_max_var_dims) + integer nx(4),ny(4) +c + i = len_trim(cfile) + if (cfile(i-2:i).eq.'.nc') then + call nchek('nf90_open-rd_missing', + & nf90_open(trim(cfile), nf90_nowrite, ncFID)) +c + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cname), ncVID)) +c + call nchek('nf90_get_att(misval)', + & nf90_get_att(ncFID, ncVID, "missing_value", misval)) +c + call nchek("nf90_close", + & nf90_close(ncFID)) + else + i = len_trim(cfile) + write(6,'(3a)') "nf90_open(",cfile(1:i-5),"," + call nchek('nf90_open-rd_missing-0000', + & nf90_open(cfile(1:i-5), nf90_nowrite, ncFID)) +c + call nchek('nf90_inq_varid', + & nf90_inq_varid(ncFID, trim(cname), ncVID)) +c + call nchek('nf90_get_att(misval)', + & nf90_get_att(ncFID, ncVID, "missing_value", misval)) +c + call nchek("nf90_close", + & nf90_close(ncFID)) + endif + return + end + + subroutine rd_bathy(n,m,h) + use mod_mom6 ! HYCOM mom6 array interface + use mod_za ! HYCOM array I/O interface + implicit none +c + integer n,m + real h(n,m) +c +c subroutine to read file for horizontal grid (depths only). +c n,m = total horizontal grid dimensions. +c h = depths (+) downward (HYCOM convention) +c + character cline*240 + character preambl(5)*79 + real hmina,hmaxa,hminb,hmaxb + integer i,j,ios +c + open (unit=9,file='regional.depth.b', + & form='formatted',status='old',action='read') + read (9, '(a79)') preambl + write(lp,'(a79)') preambl + read (9, '(a)') cline + write(lp,'(a)') trim(cline) + write(lp,'(a)') " " + call flush(lp) + i = index(cline,'=') + read (cline(i+1:),*) hminb,hmaxb + close(unit=9) +c + call zaiopf('regional.depth.a','old', 9) + call zaiord(h,ip,.false., hmina,hmaxa, 9) + call zaiocl(9) + if (abs(hmina-hminb).gt.max(abs(hmina), + & abs(hminb))*1.e-4 .or. + & abs(hmaxa-hmaxb).gt.max(abs(hmaxa), + & abs(hmaxb))*1.e-4 ) then + write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') + & 'error - .a and .b files not consistent:', + & '.a,.b min = ',hmina,hminb,hmina-hminb, + & '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb + call flush(lp) + stop + endif + do j= 1,m + do i= 1,n + if (h(i,j).gt.2.0**99) then + h(i,j) = 0.0 + endif + enddo + enddo + end + + subroutine rd_steric(n,m,ssh_mn,den_mn, flnm) + use mod_mom6 ! HYCOM mom6 array interface + use mod_za ! HYCOM array I/O interface + implicit none +c + integer n,m + real ssh_mn(n,m),den_mn(n,m) + character*256 flnm +c +c subroutine to read file for steric SSH +c n,m = total horizontal grid dimensions. +c ssh_mn = mean (steric) SSH (m) +c den_mn = mean in-situ density (kg/m^3) +c flnm = filename +c +c --- file order is consistent with HYCOM relax_ssh file +c --- den_mn (for HYCOM this would be sigma2 - thbase) +c --- ssh_mn +c --- depth - not input +c + character cline*240 + character preambl(5)*79 + real hmina,hmaxa,hminb,hmaxb + integer i,j,ios +c + open (unit=9,file=flnm(1:len_trim(flnm)-2)//'.b', + & form='formatted',status='old',action='read') + call zaiopf(flnm,'old', 9) +c +c --- den_mn +c + read (9, '(a)') cline + write(lp,'(a)') trim(cline) + call flush(lp) + i = index(cline,':') + read (cline(i+1:),*) hminb,hmaxb + call zaiord(den_mn,ip,.false., hmina,hmaxa, 9) + if (abs(hmina-hminb).gt.max(abs(hmina), + & abs(hminb))*1.e-4 .or. + & abs(hmaxa-hmaxb).gt.max(abs(hmaxa), + & abs(hmaxb))*1.e-4 ) then + write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') + & 'error - .a and .b files not consistent:', + & '.a,.b min = ',hmina,hminb,hmina-hminb, + & '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb + call flush(lp) + stop + elseif (hminb.lt.950.0 .or. hmaxb.gt.1100.0) then !sanity check + write(lp,'(/ a / a,1p2e14.6 /)') + & 'error - range not consistent with in-situ density:', + & 'min,max = ',hminb,hmaxb + call flush(lp) + stop + endif !error +c +c --- ssh_mn +c + read (9, '(a)') cline + write(lp,'(a)') trim(cline) + write(lp,'(a)') " " + call flush(lp) + i = index(cline,':') + read (cline(i+1:),*) hminb,hmaxb + call zaiord(ssh_mn,ip,.false., hmina,hmaxa, 9) + if (abs(hmina-hminb).gt.max(abs(hmina), + & abs(hminb))*1.e-4 .or. + & abs(hmaxa-hmaxb).gt.max(abs(hmaxa), + & abs(hmaxb))*1.e-4 ) then + write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') + & 'error - .a and .b files not consistent:', + & '.a,.b min = ',hmina,hminb,hmina-hminb, + & '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb + call flush(lp) + stop + elseif (hminb.lt.-5.0 .or. hmaxb.gt.5.0) then !sanity check + write(lp,'(/ a / a,1p2e14.6 /)') + & 'error - range not consistent with SSH:', + & 'min,max = ',hminb,hmaxb + endif !error +c + close(unit=9) + call zaiocl(9) + end + + subroutine rd_scp2(n,m,scp2,plat,work) + use mod_mom6 ! HYCOM mom6 array interface + use mod_za ! HYCOM array I/O interface + implicit none +c + integer n,m + real scp2(n,m),plat(n,m),work(n,m) +c +c subroutine to read file for cell size and for latitude +c n,m = total horizontal grid dimensions. +c scp2 = pscx*pscy (m) +c plat = latitude (degN) +c +c work is workspace, changed on exit +c + character cline*240 + character preambl(5)*79 + real hmina,hmaxa,hminb,hmaxb + integer i,j,ios +c + open (unit=9,file='regional.grid.b', + & form='formatted',status='old',action='read') + call zaiopf('regional.grid.a','old', 9) +c +c --- plat +c + read (9, '(a)') cline + read (9, '(a)') cline + read (9, '(a)') cline +c + read (9, '(a)') cline !plon + call zaiosk(9) +c + read (9, '(a)') cline + write(lp,*) + write(lp,'(a)') trim(cline) + call flush(lp) + i = index(cline,'=') + read (cline(i+1:),*) hminb,hmaxb +c + call zaiord(plat,ip,.false., hmina,hmaxa, 9) + if (abs(hmina-hminb).gt.max(abs(hmina), + & abs(hminb))*1.e-4 .or. + & abs(hmaxa-hmaxb).gt.max(abs(hmaxa), + & abs(hmaxb))*1.e-4 ) then + write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') + & 'error - (plat) .a and .b files not consistent:', + & '.a,.b min = ',hmina,hminb,hmina-hminb, + & '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb + call flush(lp) + stop + endif +c +c --- scpx +c + do i= 1,7 + read (9, '(a)') cline + call zaiosk(9) + enddo +c + read (9, '(a)') cline + write(lp,*) + write(lp,'(a)') trim(cline) + call flush(lp) + i = index(cline,'=') + read (cline(i+1:),*) hminb,hmaxb +c + call zaiord(work,ip,.false., hmina,hmaxa, 9) + if (abs(hmina-hminb).gt.max(abs(hmina), + & abs(hminb))*1.e-4 .or. + & abs(hmaxa-hmaxb).gt.max(abs(hmaxa), + & abs(hmaxb))*1.e-4 ) then + write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') + & 'error - (scpx) .a and .b files not consistent:', + & '.a,.b min = ',hmina,hminb,hmina-hminb, + & '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb + call flush(lp) + stop + endif +c +c --- scpy +c + read (9, '(a)') cline + write(lp,'(a)') trim(cline) + write(lp,*) + call flush(lp) + i = index(cline,'=') + read (cline(i+1:),*) hminb,hmaxb +c + call zaiord(scp2,ip,.false., hmina,hmaxa, 9) + if (abs(hmina-hminb).gt.max(abs(hmina), + & abs(hminb))*1.e-4 .or. + & abs(hmaxa-hmaxb).gt.max(abs(hmaxa), + & abs(hmaxb))*1.e-4 ) then + write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') + & 'error - (scpy) .a and .b files not consistent:', + & '.a,.b min = ',hmina,hminb,hmina-hminb, + & '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb + call flush(lp) + stop + endif +c + close(unit=9) + call zaiocl(9) +c + do j= 1,m + do i= 1,n + scp2(i,j) = scp2(i,j) * work(i,j) + enddo + enddo + end + + subroutine rd_vgrid(l,zi, cfile) + use netcdf ! NetCDF fortran 90 interface + implicit none + integer l + real zi(l+1) + character*(*) cfile +c +c subroutine to read file for fixed interface depths +c l = number of layers. +c + integer ncFID,ncVID +c + call nchek('nf90_open-rd_vgrid', + & nf90_open(trim(cfile), nf90_nowrite, ncFID)) + call nchek('nf90_inq_varid-zt_edges_ocean', + & nf90_inq_varid(ncFID,'zt_edges_ocean', ncVID)) + call nchek('nf90_get_var-zt_edges_ocean', + & nf90_get_var( ncFID, ncVID, + & zi(:))) + call nchek("nf90_close", + & nf90_close(ncFID)) +c + return + end + + subroutine nchek(cnf90,status) + use netcdf ! NetCDF fortran 90 interface + implicit none +c + character*(*), intent(in) :: cnf90 + integer, intent(in) :: status +c +c subroutine to handle NetCDF errors +c +* if (.TRUE. ) then !debug + if (.FALSE.) then !nodebug + write(6,'(a)') trim(cnf90) + endif + + if (status /= nf90_noerr) then + write(6,'(/a)') 'error from NetCDF library' + write(6,'(a/)') trim(cnf90) + write(6,'(a/)') trim(nf90_strerror(status)) + stop + end if + end subroutine nchek diff --git a/sorc/hycom_tools_for_mom6.fd/load_modules.sh b/sorc/hycom_tools_for_mom6.fd/load_modules.sh new file mode 100755 index 00000000..76645f9a --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/load_modules.sh @@ -0,0 +1,51 @@ +#!/bin/bash + +# What does this do? +# 1. Get versions of modules listed at a path. +# 2. Load them into the environment. +# + +if [[ $# -lt 1 ]]; then + echo " " + echo "Usage: " + echo "$0" "Path to RTOFS_GLO directory" + echo " " + echo "Example Inputs: " + echo " /u/santha.akella/my_ptmp/test1/RTOFS_GLO/" + echo " " + echo " " + exit 1 +fi +echo " " + +set -eu + +RTOFS_GLO_path=${1} + +source ${RTOFS_GLO_path}/versions/build.ver + +set +x +module purge +module load envvar/${envvar_ver} +module load intel/${intel_ver} +module load PrgEnv-intel/${PrgEnv_intel_ver} +module load craype/${craype_ver} +module load cray-mpich/${cray_mpich_ver} +module load cray-libsci/${cray_libsci_ver} + +module load bacio/${bacio_ver} +module load bufr/${bufr_ver} +module load g2/${g2_ver} +module load sigio/${sigio_ver} +module load sp/${sp_ver} +module load w3nco/${w3nco_ver} + +module load jasper/${jasper_ver} +module load libpng/${libpng_ver} +module load zlib/${zlib_ver} +module load netcdf/${netcdf4_ver} +module load hdf5/${hdf5_ver} + +module list +set -x + diff --git a/sorc/hycom_tools_for_mom6.fd/mod_mom6.F b/sorc/hycom_tools_for_mom6.fd/mod_mom6.F new file mode 100644 index 00000000..9479502e --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/mod_mom6.F @@ -0,0 +1,217 @@ + module mod_mom6 + implicit none +c +c --- HYCOM from mom6: array allocation interface. +c +c --- Heavily based on hycom/ALL/plot/src/mod_plot.F. +c +c --- nto = 1st dimension of mom6 t-grid +c --- mto = 2nd dimension of mom6 t-grid +c --- mvo = 2nd dimension of mom6 v-grid +c --- ii = 1st dimension of hycom (=idm) +c --- jj = 2nd dimension of hycom (=jdm) +c --- kk = actual number of layers +c --- ntracr = number of tracers +c + integer, save :: nto,ii,ii1,ii2,mto,mvo,jj,jj1,jj2,kk + integer, save :: ntracr +c +c --- input file names +c + character, save :: dpthfil*64 +c +c --- archive header +c + character, save :: ctitle(4)*80 + integer, save :: nstep,sigver +c +c --- mom6 arrays: +c + real, save, allocatable, dimension (:,:,:) :: + & s_nc,v_nc +c + real, save, allocatable, dimension (:,:) :: + & f_nc + double precision, save, allocatable, dimension (:,:) :: + & f8_nc +c + real, save, allocatable, dimension (:) :: + & zw_nc,sw_nc +c +c --- hycom arrays: +c + real, save, allocatable, dimension (:,:,:,:) :: + & trcr +c + real, save, allocatable, dimension (:,:,:) :: + & u,v,temp,saln,th3d, dp, ke +c + real, save, allocatable, dimension (:,:) :: + & ubaro,vbaro, montg,srfht,steric,dpbl,dpmixl, + & tmix,smix,thmix,umix,vmix, + & surflx,salflx,pmne, + & depths, + & covice,thkice,temice, + & kebaro,kemix, + & ttrend,strend,emnp, !never allocated + & field +c + real, save, allocatable, dimension (:) :: + & theta +c + integer, save, allocatable, dimension (:,:) :: + & ip,iq,iu,iv, + & ip0 +c +c --- module subroutines +c + contains + + subroutine mom6_alloc + implicit none +c +c --- initialize allocatable arrays. +c +c mom6 arrays +c + allocate( s_nc(nto,mto,kk+1) ) !+1 for interface tracers + allocate( v_nc(nto,mvo,kk) ) +c + allocate( f_nc(nto,mto) ) +c + allocate( zw_nc(kk+1) ) + allocate( sw_nc(kk+1) ) +c +c hycom arrays. +c + ii1 = ii - 1 + ii2 = ii - 2 + jj1 = jj - 1 + jj2 = jj - 2 +c + allocate( u(ii,jj,kk) ) + allocate( v(ii,jj,kk) ) + allocate( temp(ii,jj,kk) ) + allocate( saln(ii,jj,kk) ) + allocate( th3d(ii,jj,kk) ) + allocate( dp(ii,jj,kk) ) + allocate( ke(ii,jj,kk) ) !may not be needed +c + if (ntracr.gt.0) then + allocate( trcr(ii,jj,kk,ntracr) ) + endif +c + allocate( ubaro(ii,jj) ) + allocate( vbaro(ii,jj) ) + allocate( montg(ii,jj) ) + allocate( srfht(ii,jj) ) + allocate( steric(ii,jj) ) + allocate( dpbl(ii,jj) ) + allocate( dpmixl(ii,jj) ) + allocate( tmix(ii,jj) ) + allocate( smix(ii,jj) ) + allocate( thmix(ii,jj) ) + allocate( umix(ii,jj) ) + allocate( vmix(ii,jj) ) + allocate( pmne(ii,jj) ) + allocate( salflx(ii,jj) ) + allocate( surflx(ii,jj) ) + allocate( covice(ii,jj) ) + allocate( thkice(ii,jj) ) + allocate( temice(ii,jj) ) + allocate( kebaro(ii,jj) ) + allocate( kemix(ii,jj) ) +c + allocate( depths(0:ii,0:jj) ) + allocate( ip0(0:ii,0:jj) ) +c + allocate( ip(ii,jj) ) + allocate( iq(ii,jj) ) + allocate( iu(ii,jj) ) + allocate( iv(ii,jj) ) +c + allocate( theta(kk) ) + + end subroutine mom6_alloc + + subroutine mom6_alloc_tide + implicit none +c +c --- initialize allocatable arrays for ssh. +c +c mom6 arrays +c + allocate( f_nc(nto,mvo) ) !mvo might be mto+1 +c +c hycom arrays. +c + ii1 = ii - 1 + ii2 = ii - 2 + jj1 = jj - 1 + jj2 = jj - 2 +c + allocate( montg(ii,jj) ) + allocate( srfht(ii,jj) ) + allocate( steric(ii,jj) ) + allocate( ubaro(ii,jj) ) + allocate( vbaro(ii,jj) ) +c + allocate( depths(0:ii,0:jj) ) + allocate( ip0(0:ii,0:jj) ) +c + allocate( ip(ii,jj) ) + allocate( iu(ii,jj) ) + allocate( iv(ii,jj) ) + allocate( iq(ii,jj) ) + + end subroutine mom6_alloc_tide + + subroutine mom6_alloc_field + implicit none +c +c --- initialize allocatable arrays for 2field. +c +c mom6 arrays +c + allocate( f_nc(nto,mto) ) +c +c hycom arrays. +c + allocate( ip(ii,jj) ) + allocate( field(ii,jj) ) + + end subroutine mom6_alloc_field + + subroutine mom6_alloc_field_3d + implicit none +c +c --- initialize allocatable arrays for 3d field. +c +c mom6 arrays +c + allocate( s_nc(nto,mto,kk) ) +c +c hycom arrays. +c + allocate( ip(ii,jj) ) + allocate( field(ii,jj) ) + + end subroutine mom6_alloc_field_3d + + subroutine mom6_alloc_field8 + implicit none +c +c --- initialize allocatable arrays for 2field. +c +c mom6 arrays +c + allocate( f8_nc(nto,mto) ) +c +c hycom arrays. +c + allocate( ip(ii,jj) ) + allocate( field(ii,jj) ) + + end subroutine mom6_alloc_field8 + + end module mod_mom6 diff --git a/sorc/hycom_tools_for_mom6.fd/mod_ncio.f90 b/sorc/hycom_tools_for_mom6.fd/mod_ncio.f90 new file mode 100644 index 00000000..c2e51d15 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/mod_ncio.f90 @@ -0,0 +1,595 @@ +Module mod_ncio +use netcdf + +!---------------------------------------- +! extracted from TSIS package +! version 0.1 -- 2013/03/02 +! Ashwanth Srinivasan +! Tendral LLC, Key Biscayne +! Modified by Alexandra Bozec (03/2021) +! for HYCOM-tools/mom6, FSU/COAPS +! --------------------------------------- +! a minimal module for handling netcdf io. Implements the most +! commonly used functionalities as a set of subroutines wraping all lower +! level netcdf calls. All call take in a file name and a file id returned by +! the netcdf lib. +! nciocf - creates a netcdf file +! nciopn - opens a netcdf file +! nciocl - closes a netcdf file +! ncioed - end definition mode and enters data mode +! nciocv - creates a variable +! nciorv - reads a varaible - interface to nciordi1-4,nciordr1-4,nciordd1-4 +! nciowv - writes a varaible - interface to nciowri1-4,nciowrr1-4,nciowrd1-4 +! ncioin - inquires about a variable +! nciowa - put variable attributes +! nciora - get variable attributes +! nciock - error handling routine + +implicit none +! define kinds +integer, parameter :: cl = 120, & + i4 = selected_int_kind(9), & + i8 = selected_int_kind(18),& + bl = kind(.true.), & + r4 = selected_real_kind(6),& + r8 = selected_real_kind(12) + +interface nciorv + module procedure nciorvi1 + module procedure nciorvr1 + module procedure nciorvd1 + module procedure nciorvi2 + module procedure nciorvr2 + module procedure nciorvI3 + module procedure nciorvr3 + module procedure nciorvr4 + module procedure nciorvca +end interface + + +interface nciowv + module procedure nciowvi1 + module procedure nciowvr1 +! module procedure nciowvd1 + module procedure nciowvI2 + module procedure nciowvr2 + module procedure nciowvI3 + module procedure nciowvr3 + module procedure nciowvI4 + module procedure nciowvr4 +end interface + +interface nciowa + module procedure nciowac + module procedure nciowar + module procedure nciowgac + module procedure nciowgai + module procedure nciowgar +end interface + +interface nciora + module procedure nciordgai + module procedure nciordgar + module procedure nciordar +end interface + + +contains + +!### +subroutine nciocf(filename,fid) +implicit none +character(len=*), intent(in) :: filename +integer(i4), intent(out) :: fid +integer(i4) :: ncid +!!Alex call nciock(filename,nf90_create(path = trim(filename),cmode=NF90_64BIT_OFFSET, ncid = fid)) +call nciock(filename,nf90_create(path = trim(filename),cmode=or(nf90_clobber, & + or(nf90_hdf5, & + nf90_classic_model)), ncid = fid)) + +end subroutine + +!### +subroutine ncioed(filename,fid) +implicit none +character(len=*), intent(in) :: filename +integer(i4), intent(in) :: fid +call nciock(filename,nf90_enddef(fid)) +end subroutine + +!### +subroutine nciorc(filename,fid) + implicit none + character(len=*), intent(in) :: filename + integer(i4), intent(in) :: fid + call nciock(filename,nf90_redef(fid)) +end subroutine + + + +!### +subroutine nciodd(filename,fid,dimname,dimlen) +implicit none +character(len=*), intent(in) :: filename +integer(i4), intent(in) :: fid +character(len=*), intent(in) :: dimname +integer(i4), intent(in) :: dimlen +integer(i4) :: dimid +call nciock(filename,nf90_def_dim(fid,trim(dimname), dimlen,dimid)) +end subroutine + +!### +subroutine nciocv(filename,fid,vname,vtype,vdims) +implicit none +character(len=*), intent(in) :: filename +integer(i4), intent(in) :: fid +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: vtype +integer(i4), intent(in) :: vdims(:) +integer(i4) :: varid +call nciock(filename,nf90_def_var(fid,trim(vname),vtype,vdims, varid)) +end subroutine + +!### +subroutine nciopn(filename,fid) +implicit none +character(len=*), intent(in) :: filename +integer(i4), intent(out) :: fid +call nciock(filename,nf90_open(trim(filename),nf90_Write,fid)) +end subroutine nciopn + +!### +subroutine nciocl(filename, fid) +implicit none +character(len=*), intent(in) :: filename +integer(i4), intent(in) :: fid +call nciock(filename,nf90_close(fid)) +end subroutine nciocl + +!### +subroutine ncioin(filename,fid,vname,nval) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +integer(i4), intent(out) :: nval +integer(i4) :: vardimid +call nciock(filename,nf90_inq_dimid(fid,trim(vname), vardimid)) +call nciock(filename,nf90_inquire_dimension(fid,vardimid, len = nval)) +end subroutine + + +SUBROUTINE ncioiv(ncfname,ncId,varName,varid) +CHARACTER (LEN=*), INTENT(IN) :: ncfname +CHARACTER (LEN=*), INTENT(IN) :: varName +INTEGER(i4), INTENT(IN) :: ncId +INTEGER(i4) :: varID,status + +status=nf90_inq_varid(ncid,trim(varName),varid) +varid=status +end subroutine + +!### +subroutine nciorvca(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +character(len=*), intent(inout) :: values +integer(i4),intent(in),optional :: start(2) +integer(i4),intent(in),optional :: count(2) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + + +!1d integer variable +subroutine nciorvi1(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +integer(i4),dimension(:), intent(inout) :: values +integer(i4),intent(in),optional :: start(1) +integer(i4),intent(in),optional :: count(1) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + +!### +!1d real variable +subroutine nciorvr1(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4),dimension(:), intent(inout) :: values +integer(i4),intent(in),optional :: start(1) +integer(i4),intent(in),optional :: count(1) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + +!### +! 1d double variable +subroutine nciorvd1(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r8), dimension(:), intent(inout) :: values +integer(i4), intent(in),optional :: start(1) +integer(i4), intent(in),optional :: count(1) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + +!### + +!2d int variable +subroutine nciorvI2(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +integer(i4), dimension(:,:), intent(inout) :: values +integer(i4), intent(in),optional :: start(2) +integer(i4), intent(in),optional :: count(2) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + +!2d real variable +subroutine nciorvr2(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4), dimension(:,:), intent(inout) :: values +integer(i4), intent(in),optional :: start(2) +integer(i4), intent(in),optional :: count(2) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + +!### 3d int variable +subroutine nciorvI3(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +integer(r4), dimension(:,:,:), intent(inout) :: values +integer(i4), intent(in),optional :: start(3) +integer(i4), intent(in),optional :: count(3) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + +!### 3d real variable +subroutine nciorvr3(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4), dimension(:,:,:), intent(inout) :: values +integer(i4), intent(in),optional :: start(3) +integer(i4), intent(in),optional :: count(3) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_get_var(fid, vid, values)) +endif +end subroutine + +!### 4d real variable +subroutine nciorvr4(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4), dimension(:,:,:,:), intent(inout) :: values +integer(i4), intent(in),optional :: start(4) +integer(i4), intent(in),optional :: count(4) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_get_var(fid, vid, values,start,count)) + else + call nciock(filename,nf90_get_var(fid, vid, values)) + endif + end subroutine + +!### write a 1d integer variable + +!### write a 1d real variable +subroutine nciowvi1(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +integer(i4),dimension(:), intent(in) :: values +integer(i4),intent(in),optional :: start(1) +integer(i4),intent(in),optional :: count(1) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_put_var(fid, vid, values)) +endif +end subroutine + + +!### write a 1d real variable +subroutine nciowvr1(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4),dimension(:), intent(in) :: values +integer(i4),intent(in),optional :: start(1) +integer(i4),intent(in),optional :: count(1) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_put_var(fid, vid, values)) +endif +end subroutine + +!### write 2d int variable +subroutine nciowvi2(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +integer(i4), dimension(:,:), intent(in) :: values +integer(i4), intent(in),optional :: start(2) +integer(i4), intent(in),optional :: count(2) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_put_var(fid, vid, values)) +endif +end subroutine +!### write 2d real variable +subroutine nciowvr2(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4), dimension(:,:), intent(in) :: values +integer(i4), intent(in),optional :: start(2) +integer(i4), intent(in),optional :: count(2) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_put_var(fid, vid, values)) +endif +end subroutine + +!### write 3d int variable +subroutine nciowvI3(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +integer(r4), dimension(:,:,:), intent(in) :: values +integer(i4), intent(in),optional :: start(3) +integer(i4), intent(in),optional :: count(3) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_put_var(fid, vid, values)) +endif +end subroutine + +!### write a 3d real variable +subroutine nciowvr3(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4), dimension(:,:,:), intent(in) :: values +integer(i4), intent(in),optional :: start(3) +integer(i4), intent(in),optional :: count(3) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid, values,start,count)) +else + call nciock(filename,nf90_put_var(fid, vid, values)) +endif +end subroutine + +!### write a 4d int variable +subroutine nciowvI4(filename,fid,vname,values,start,count) + implicit none + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: vname + integer(i4), intent(in) :: fid + integer(r4), dimension(:,:,:,:), intent(in) :: values + integer(i4), intent(in),optional :: start(4) + integer(i4), intent(in),optional :: count(4) + integer(i4) :: vid + call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) + if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid,values,start,count)) + else + call nciock(filename,nf90_put_var(fid, vid, values)) + endif +end subroutine + +!### write a 4d real variable +subroutine nciowvr4(filename,fid,vname,values,start,count) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +integer(i4), intent(in) :: fid +real(r4), dimension(:,:,:,:), intent(in) :: values +integer(i4), intent(in),optional :: start(4) +integer(i4), intent(in),optional :: count(4) +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +if (present(start)) then + call nciock(filename,nf90_put_var(fid, vid, values,start,count)) + else + call nciock(filename,nf90_put_var(fid, vid, values)) + endif +end subroutine + +!### write/put variable attributes +subroutine nciowac(filename,fid,vname,name,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +character(len=*), intent(in) :: name +character(len=*), intent(in) :: value +integer(i4), intent(in) :: fid +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +call nciock(filename,nf90_put_att(fid, vid, name,value)) +end subroutine + +!### +subroutine nciowar(filename,fid,vname,name,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: vname +character(len=*), intent(in) :: name +real(r4), intent(in) :: value +integer(i4), intent(in) :: fid +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid,trim(vname), vid)) +call nciock(filename,nf90_put_att(fid, vid, name,value)) +end subroutine + +!### +subroutine nciowgar(filename,fid,name,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: name +real(r4), intent(in) :: value +integer(i4), intent(in) :: fid +call nciock(filename,nf90_put_att(fid, NF90_GLOBAL, name,value)) +end subroutine + +!### +subroutine nciowgac(filename,fid,name,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: name +character(len=*), intent(in) :: value +integer(i4), intent(in) :: fid +call nciock(filename,nf90_put_att(fid, NF90_GLOBAL, name,value)) +end subroutine + +!### +subroutine nciowgai(filename,fid,name,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: name +integer(i4), intent(in) :: value +integer(i4), intent(in) :: fid +call nciock(filename,nf90_put_att(fid, NF90_GLOBAL, name,value)) +end subroutine + +!### +!### +subroutine nciordgai(filename,fid,name,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: name +integer(i4), intent(inout) :: value +integer(i4), intent(in) :: fid +call nciock(filename,nf90_get_att(fid, NF90_GLOBAL, name,value)) +end subroutine + +!### +subroutine nciordgar(filename,fid,name,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: name +REAL(r4), intent(inout) :: value +integer(i4), intent(in) :: fid +call nciock(filename,nf90_get_att(fid, NF90_GLOBAL, name,value)) +end subroutine + +!### +subroutine nciordar(filename,fid,namevar,nameatt,value) +implicit none +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: namevar +character(len=*), intent(in) :: nameatt +REAL(r4), intent(inout) :: value +integer(i4), intent(in) :: fid +integer(i4) :: vid +call nciock(filename,nf90_inq_varid(fid, namevar, vid)) +call nciock(filename,nf90_get_att(fid, vid, nameatt, value)) +end subroutine + +!### +subroutine nciock(filename,istat) +use netcdf +implicit none + +!This routine provides a simple interface to netCDF error message routine. +!provides a filename and an error message + +character(len=*), intent(in) :: filename +integer,intent(in) :: istat + +if (istat /= nf90_noerr) then + write (*,*) 'Error operating on netCDF file: ', adjustl(trim(filename)) + write (*,*) trim(nf90_strerror(istat)) + stop +endif + +end subroutine nciock + + +end module mod_ncio diff --git a/sorc/hycom_tools_for_mom6.fd/mod_xc.F b/sorc/hycom_tools_for_mom6.fd/mod_xc.F new file mode 100644 index 00000000..c031bf42 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/mod_xc.F @@ -0,0 +1,323 @@ + module mod_xc + implicit none +c +c --- HYCOM communication interface. +c --- A subset of the serial interface for setup only. +c +c --- tital array dimensions + integer, public, save :: idm,jdm +c +c --- halo size always zero for setup + integer nbdy + parameter (nbdy=0) +c +c --- line printer unit (stdout) + integer lp + common/linepr/ lp + save /linepr/ +c +c --- tile number (counting from 1) + integer, public, save :: mnproc +c +c --- xcsync stdout flushing options + logical, public, parameter :: flush_lp=.true., + & no_flush=.false. +c +c --- private timer variables, see xctmri + character*6, private, dimension(97), save :: cc + integer, private, dimension(97), save :: nc + real*8, private, dimension(97), save :: tc,t0 +c +c --- actual module subroutines + contains + + subroutine xcspmd + implicit none +c +c********** +c* +c 1) initialize data structures that identify the domain and tiles. +c +c 2) data structures: +c idm - 1st total array dimension +c jdm - 2nd total array dimension +c mnproc - 1-D node index +c +c 3) Total array dimensions from regional.grid.b +c* +c********** +c + character cvarin*6 +c +c shared memory version, mnproc=1. +c + mnproc = 1 + lp = 6 +c +c total array dimensions from regional.grid.b +c + open(unit=11,file='regional.grid.b',form='formatted', + & status='old',action='read') +c + read( 11,*) idm,cvarin + if (cvarin.ne.'idm ') then + write(lp,*) + write(lp,*) 'error in xcspmd - regional.grid.b input ',cvarin, + & ' but should be idm ' + write(lp,*) + stop + endif + read( 11,*) jdm,cvarin + if (cvarin.ne.'jdm ') then + write(lp,*) + write(lp,*) 'error in xcspmd - regional.grid.b input ',cvarin, + & ' but should be jdm ' + write(lp,*) + call flush(lp) + stop + endif +c + write(lp,'(/ a,2i5 /)') 'xcspmd: idm,jdm =',idm,jdm +c + close(unit=11) +c +c initialize timers. +c + call xctmri + return + end subroutine xcspmd + + subroutine xcstop(cerror) + implicit none +c + character*(*), intent(in) :: cerror +c +c********** +c* +c 1) stop all processes. +c +c 2) all processes must call this routine. +c use 'xchalt' for emergency stops. +c +c 3) parameters: +c name type usage description +c ---------- ---------- ------- ---------------------------- +c cerror char*(*) input error message +c* +c********** +c +c print active timers. +c + call xctmrp +c +c shared memory version, just stop. +c + if (cerror.ne.' ') then + write(lp,*) '**************************************************' + write(lp,*) cerror + write(lp,*) '**************************************************' + call flush(lp) + endif + stop '(xcstop)' + end subroutine xcstop + + subroutine xcsync(lflush) + implicit none +c + logical, intent(in) :: lflush +c +c********** +c* +c 1) barrier, no processor exits until all arrive (and flush stdout). +c +c 2) some MPI implementations only flush stdout as a collective +c operation, and hence the lflush=.true. option to flush stdout. +c +c 3) Only one processor, so the barrier is a no-op in this case. +c* +c********** +c + if (lflush) then + call flush(lp) + endif + return + end subroutine xcsync + + subroutine xctmri + implicit none +c +c +c********** +c* +c 1) initialize timers. +c +c 2) timers 1:32 are for message passing routines, +c timers 33:80 are for general hycom routines, +c timers 81:96 are for user selected routines. +c timer 97 is the total time. +c +c 3) call xctmri to initialize timers (called in xcspmd), +c call xctmr0(n) to start timer n, +c call xctmr1(n) to stop timer n and add event to timer sum, +c call xctnrn(n,cname) to register a name for timer n, +c call xctmrp to printout timer statistics (called by xcstop). +c* +c********** +c + integer i +c + real*8 zero8 + parameter (zero8=0.0) +c + do 110 i= 1,97 + cc(i) = ' ' + nc(i) = 0 + tc(i) = zero8 + 110 continue +c + call xctmrn(97,'total ') + call xctmr0(97) + return + end subroutine xctmri + + subroutine xctmr0(n) + implicit none +c + integer, intent(in) :: n +c +c********** +c* +c 1) start timer n. +c +c 2) parameters: +c name type usage description +c ---------- ---------- ------- ---------------------------- +c n integer input timer number +c* +c********** +c + real*8 wtime +c +#if defined(DEBUG_TIMER) + if (n.gt.24 .and. cc(n).ne.' ') then + write(lp,*) 'call ',cc(n) + call flush(lp) + endif +#endif + t0(n) = wtime() + return + end subroutine xctmr0 + + subroutine xctmr1(n) + implicit none +c + integer, intent(in) :: n +c +c********** +c* +c 1) add time since call to xctim0 to timer n. +c +c 2) parameters: +c name type usage description +c ---------- ---------- ------- ---------------------------- +c n integer input timer number +c* +c********** +c + real*8 wtime +c + nc(n) = nc(n) + 1 + tc(n) = tc(n) + (wtime() - t0(n)) +#if defined(DEBUG_TIMER) + if (n.gt.24 .and. cc(n).ne.' ') then + write(lp,*) 'exit ',cc(n) + call flush(lp) + endif +#endif + return + end subroutine xctmr1 + + subroutine xctmrn(n,cname) + implicit none +c + character*6, intent(in) :: cname + integer, intent(in) :: n +c +c********** +c* +c 1) register name of timer n. +c +c 2) parameters: +c name type usage description +c ---------- ---------- ------- ---------------------------- +c n integer input timer number +c cname char*(8) input timer name +c* +c********** +c + cc(n) = cname + return + end subroutine xctmrn + + subroutine xctmrp + implicit none +c +c********** +c* +c 1) print all active timers. +c +c 2) on exit all timers are reset to zero. +c* +c********** +c + integer i +c + real*8 zero8 + parameter (zero8=0.0) +c +c get total time. +c + call xctmr1(97) +c +c print timers. +c + write(lp,6000) + do i= 1,97 + if (nc(i).ne.0) then + if (cc(i).ne.' ') then + write(lp,6100) cc(i),nc(i),tc(i),tc(i)/nc(i) + else + write(lp,6150) i, nc(i),tc(i),tc(i)/nc(i) + endif + endif + enddo + write(lp,6200) + call flush(lp) +c +c reset timers to zero. +c + do i= 1,97 + nc(i) = 0 + tc(i) = zero8 + enddo +c +c start a new total time measurement. +c + call xctmr0(97) + return +c + 6000 format(/ / + + 4x,' timer statistics ' / + + 4x,'------------------' /) + 6100 format(5x,a6, + + ' calls =',i9, + + ' time =',f11.5, + + ' time/call =',f14.8) + 6150 format(5x,' #',i2, + + ' calls =',i9, + + ' time =',f11.5, + + ' time/call =',f14.8) + 6200 format(/ /) + end subroutine xctmrp + + end module mod_xc diff --git a/sorc/hycom_tools_for_mom6.fd/mod_za.F b/sorc/hycom_tools_for_mom6.fd/mod_za.F new file mode 100644 index 00000000..eed3a47f --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/mod_za.F @@ -0,0 +1,1026 @@ + module mod_za + use mod_xc ! HYCOM communication API +c + implicit none +c +c --- HYCOM I/O interface. +c --- Serial version, for setup only. +c + integer, save, private :: iarec(999) + real*4, save, private, allocatable :: w(:) +c +c n2drec = size of output 2-d array, multiple of 4096 +c spval = data void marker, 2^100 or about 1.2676506e30 +c + integer, save, private :: n2drec + real*4, private, parameter :: spval=2.0**100 +c + private zaiordd,zaiowrd + + contains + +c +c----------------------------------------------------------------------- +c +c machine dependent I/O routines. +c single processor version, contained in mod_za. +c +c author: Alan J. Wallcraft, NRL. +c +c----------------------------------------------------------------------- +c + subroutine zaiopn(cstat, iaunit) + implicit none +c + integer, intent(in) :: iaunit + character*(*), intent(in) :: cstat +c +c********** +c* +c 1) machine specific routine for opening a file for array i/o. +c +c must call zaiost before first call to zaiopn. +c see also 'zaiope' and 'zaiopf'. +c +c 2) the filename is taken from the environment variable FORxxxA, +c where xxx = iaunit, with default fort.xxxa. +c +c array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c cstat indicates the file type, it can be 'scratch', 'old', or +c 'new'. +c all i/o to iaunit must be performed by zaiord and zaiowr. +c the file should be closed using zaiocl. +c* +c********** +c + integer ios,nrecl + character cfile*256,cenv*7 + character cact*9 +#if defined(TIMER) +c + call xctmr0(16) +#endif +c +c test file state. +c + if (iarec(iaunit).ne.-1) then + write(6,9000) iaunit + call flush(6) + stop + endif +c +c get filename. +c + write(cenv,1000) iaunit + cfile = ' ' + call getenv(cenv,cfile) + if (cfile.eq.' ') then + write(cfile,1100) iaunit + endif +* write(6,*) 'zaiopn - iaunit = ',iaunit +* call flush(6) +c +c open file. +c + inquire(iolength=nrecl) w +c + if (cstat.eq.'OLD' .or. + & cstat.eq.'old' ) then + cact = 'READ' + elseif (cstat.eq.'NEW' .or. + & cstat.eq.'new' ) then + cact = 'WRITE' + else + cact = 'READWRITE' + endif +#if defined(YMP) + if (mod(nrecl,16384).eq.0 .and. nrecl.gt.16384*4) then + call asnunit(iaunit+1000,'-F syscall -N ieee',ios) + else + call asnunit(iaunit+1000,'-F cachea:8:16:2 -N ieee',ios) + endif + if (ios.ne.0) then + write(6,9050) iaunit + write(6,*) 'ios = ',ios + call flush(6) + stop + endif +#endif + if (cstat.eq.'scratch' .or. + & cstat.eq.'SCRATCH' ) then + open(unit=iaunit+1000, + & form='unformatted', status='scratch', + & access='direct', recl=nrecl, action=cact, iostat=ios) + else + open(unit=iaunit+1000, file=cfile, + & form='unformatted', status=cstat, + & access='direct', recl=nrecl, action=cact, iostat=ios) + endif + if (ios.ne.0) then + write(6,9100) iaunit + write(6,*) 'ios = ',ios + call flush(6) + stop + endif + iarec(iaunit) = 0 +#if defined(TIMER) +c + call xctmr1(16) +#endif + return +c + 1000 format('FOR',i3.3,'A') + 1100 format('fort.',i3.3,'a') + 9000 format(/ /10x,'error in zaiopn - array I/O unit ', + & i3,' is not marked as available.'/ /) +#if defined(YMP) + 9050 format(/ /10x,'error in zaiopn - can''t asnunit ',i3, + & ', for array I/O.'/ /) +#endif + 9100 format(/ /10x,'error in zaiopn - can''t open unit ',i3, + & ', for array I/O.'/ /) + end subroutine zaiopn + + subroutine zaiope(cenv,cstat, iaunit) + implicit none +c + integer, intent(in) :: iaunit + character*(*), intent(in) :: cenv,cstat +c +c********** +c* +c 1) machine specific routine for opening a file for array i/o. +c +c must call zaiost before first call to zaiope. +c see also 'zaiopn' and 'zaiopf'. +c +c 2) the filename is taken from environment variable 'cenv'. +c +c array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c cstat indicates the file type, it can be 'scratch', 'old', or +c 'new'. +c all i/o to iaunit must be performed by zaiord and zaiowr. +c arrays passed to these routines must conform to 'h'. +c the file should be closed using zaiocl. +c* +c********** +c + integer ios,nrecl + character cfile*256 + character cact*9 +#if defined(TIMER) +c + call xctmr0(16) +#endif +c +c test file state. +c + if (iarec(iaunit).ne.-1) then + write(6,9000) iaunit + call flush(6) + stop + endif +c +c get filename. +c + cfile = ' ' + call getenv(cenv,cfile) + if (cfile.eq.' ') then + write(6,9300) cenv(1:len_trim(cenv)) + write(6,*) 'iaunit = ',iaunit + call flush(6) + stop + endif +c +c open file. +c +* write(6,*) 'zaiope - iaunit = ',iaunit +* call flush(6) +* + inquire(iolength=nrecl) w +c + if (cstat.eq.'OLD' .or. + & cstat.eq.'old' ) then + cact = 'READ' + elseif (cstat.eq.'NEW' .or. + & cstat.eq.'new' ) then + cact = 'WRITE' + else + cact = 'READWRITE' + endif +c +#if defined(YMP) + if (mod(nrecl,16384).eq.0 .and. nrecl.gt.16384*4) then + call asnunit(iaunit+1000,'-F syscall -N ieee',ios) + else + call asnunit(iaunit+1000,'-F cachea:8:16:2 -N ieee',ios) + endif + if (ios.ne.0) then + write(6,9050) iaunit,cfile(1:len_trim(cfile)) + write(6,*) 'ios = ',ios + write(6,*) 'cenv = ',cenv(1:len_trim(cenv)) + call flush(6) + stop + endif +#endif + open(unit=iaunit+1000, file=cfile, + & form='unformatted', status=cstat, + & access='direct', recl=nrecl, action=cact, iostat=ios) + if (ios.ne.0) then + write(6,9100) iaunit,cfile(1:len_trim(cfile)) + write(6,*) 'ios = ',ios + write(6,*) 'cenv = ',cenv(1:len_trim(cenv)) + call flush(6) + stop + endif + iarec(iaunit) = 0 +#if defined(TIMER) +c + call xctmr1(16) +#endif + return +c + 9000 format(/ /10x,'error in zaiope - array I/O unit ', + & i3,' is not marked as available.'/ /) +#if defined(YMP) + 9050 format(/ /10x,'error in zaiope - can''t asnunit ',i3, + & ', for array I/O.' / + & 10x,'cfile = ',a/ /) +#endif + 9100 format(/ /10x,'error in zaiope - can''t open unit ',i3, + & ', for array I/O.' / + & 10x,'cfile = ',a/ /) + 9300 format(/ /10x,'error in zaiope - environment variable ',a, + & ' not defined'/ /) + end subroutine zaiope + + subroutine zaiopf(cfile,cstat, iaunit) + implicit none +c + integer, intent(in) :: iaunit + character*(*), intent(in) :: cfile,cstat +c +c********** +c* +c 1) machine specific routine for opening a file for array i/o. +c +c must call zaiost before first call to zaiopf. +c see also 'zaiopn' and 'zaiope'. +c +c 2) the filename is taken from 'cfile'. +c +c array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c cstat indicates the file type, it can be 'scratch', 'old', or +c 'new'. +c all i/o to iaunit must be performed by zaiord and zaiowr. +c arrays passed to these routines must conform to 'h'. +c the file should be closed using zaiocl. +c* +c********** +c + integer ios,nrecl + character cact*9 +#if defined(TIMER) +c + call xctmr0(16) +#endif +c +c test file state. +c + if (iarec(iaunit).ne.-1) then + write(6,9000) iaunit + call flush(6) + stop + endif +c +c open file. +c +* write(6,*) 'zaiopf - iaunit = ',iaunit +* call flush(6) +* + inquire(iolength=nrecl) w +c + if (cstat.eq.'OLD' .or. + & cstat.eq.'old' ) then + cact = 'READ' + elseif (cstat.eq.'NEW' .or. + & cstat.eq.'new' ) then + cact = 'WRITE' + else + cact = 'READWRITE' + endif +c +#if defined(YMP) + if (mod(nrecl,16384).eq.0 .and. nrecl.gt.16384*4) then + call asnunit(iaunit+1000,'-F syscall -N ieee',ios) + else + call asnunit(iaunit+1000,'-F cachea:8:16:2 -N ieee',ios) + endif + if (ios.ne.0) then + write(6,9050) iaunit,cfile(1:len_trim(cfile)) + write(6,*) 'ios = ',ios + call flush(6) + stop + endif +#endif + open(unit=iaunit+1000, file=cfile, + & form='unformatted', status=cstat, + & access='direct', recl=nrecl, action=cact, iostat=ios) + if (ios.ne.0) then + write(6,9100) iaunit,cfile(1:len_trim(cfile)) + write(6,*) 'ios = ',ios + call flush(6) + stop + endif + iarec(iaunit) = 0 +#if defined(TIMER) +c + call xctmr1(16) +#endif + return +c + 9000 format(/ /10x,'error in zaiopf - array I/O unit ', + & i3,' is not marked as available.'/ /) +#if defined(YMP) + 9050 format(/ /10x,'error in zaiopf - can''t asnunit ',i3, + & ', for array I/O.' / + & 10x,'cfile = ',a/ /) +#endif + 9100 format(/ /10x,'error in zaiopf - can''t open unit ',i3, + & ', for array I/O.' / + & 10x,'cfile = ',a/ /) + end subroutine zaiopf + + subroutine zaiopi(lopen, iaunit) + implicit none +c + logical, intent(out) :: lopen + integer, intent(in) :: iaunit +c +c********** +c* +c 1) is an array i/o unit open? +c +c 2) must call zaiost before first call to zaiopi. +c* +c********** +c + lopen = iarec(iaunit).ne.-1 + return + end subroutine zaiopi + + subroutine zaiost + implicit none +c +c********** +c* +c 1) machine specific routine for initializing array i/o. +c +c 2) see also zaiopn, zaiord, zaiowr, and zaiocl. +c* +c********** +c +c n2drec = size of output 2-d array, multiple of 4096 +c + n2drec = ((idm*jdm+4095)/4096)*4096 +c +c initialize I/O buffer +c + allocate( w(n2drec) ) +c +c initialize record counters +c + iarec(:) = -1 +#if defined(TIMER) +c +c initialize timers. +c + call xctmrn(16,'zaio**') + call xctmrn(17,'zaiord') + call xctmrn(18,'zaiowr') +#endif + return + end subroutine zaiost + + subroutine zaiocl(iaunit) + implicit none +c + integer, intent(in) :: iaunit +c +c********** +c* +c 1) machine specific routine for array i/o file closing. +c +c must call zaiopn for this array unit before calling zaiocl. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c* +c********** +c + integer ios +#if defined(TIMER) +c + call xctmr0(16) +#endif +c +* write(6,*) 'zaiocl - iaunit = ',iaunit +* call flush(6) + if (iarec(iaunit).lt.0) then + write(6,9000) iaunit + call flush(6) + stop + endif +c + close(unit=iaunit+1000, status='keep') +#if defined(T3E) || defined(YMP) || defined(X1) + call asnunit(iaunit+1000,'-R',ios) +#endif + iarec(iaunit) = -1 +#if defined(TIMER) +c + call xctmr1(16) +#endif + return +c + 9000 format(/ /10x,'error in zaiocl - array I/O unit ', + & i3,' is not marked as open.'/ /) + end subroutine zaiocl + + subroutine zaiofl(iaunit) + implicit none +c + integer, intent(in) :: iaunit +c +c********** +c* +c 1) machine specific routine for array i/o buffer flushing. +c +c must call zaiopn for this array unit before calling zaiocl. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c* +c********** +c + integer irlen + character cfile*256 +#if defined(TIMER) +c + call xctmr0(16) +#endif +c + if (iarec(iaunit).lt.0) then + write(6,9000) iaunit + call flush(6) + stop + endif +c + inquire(unit=iaunit+1000, name=cfile, recl=irlen) + close( unit=iaunit+1000, status='keep') + open( unit=iaunit+1000, file=cfile, form='unformatted', + & access='direct', recl=irlen) +#if defined(TIMER) +c + call xctmr1(16) +#endif + return +c + 9000 format(/ /10x,'error in zaiofl - array I/O unit ', + & i3,' is not marked as open.'/ /) + end subroutine zaiofl + + subroutine zaiorw(iaunit) + implicit none +c + integer, intent(in) :: iaunit +c +c********** +c* +c 1) machine specific routine for array i/o file rewinding. +c +c must call zaiopn for this array unit before calling zaiocl. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c* +c********** +#if defined(TIMER) +c + call xctmr0(16) +#endif +c + if (iarec(iaunit).lt.0) then + write(6,9000) iaunit + call flush(6) + stop + endif +c + iarec(iaunit) = 0 +* write(6,*) 'zaiorw - iaunit,rec = ',iaunit,iarec(iaunit) +* call flush(6) +#if defined(TIMER) +c + call xctmr1(16) +#endif + return +c + 9000 format(/ /10x,'error in zaiorw - array I/O unit ', + & i3,' is not marked as open.'/ /) + end subroutine zaiorw + + subroutine zaiord3(h, l, mask,lmask, hmin,hmax, iaunit) + implicit none +c + logical, intent(in) :: lmask + integer, intent(in) :: l,iaunit + integer, dimension (1:idm,1:jdm), + & intent(in) :: mask +#if defined(REAL4) + real*4, intent(out) :: hmin(l),hmax(l) + real*4, dimension (1:idm,1:jdm,l), + & intent(out) :: h +#else + real, intent(out) :: hmin(l),hmax(l) + real, dimension (1:idm,1:jdm,l), + & intent(out) :: h +#endif +c +c********** +c* +c 1) machine specific routine for 3-d array reading. +c +c must call zaiopn for this array unit before calling zaiord. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c the array, 'h', must conform to that passed in the associated +c call to zaiopn. +c +c 4) hmin,hmax are returned as the minimum and maximum value in the +c array, ignoring array elements set to 2.0**100. +c if lmask==.true. the range is calculated only where mask.ne.0, +c with all other values unchanged in h on exit. It is then an +c error if mask.ne.0 anywhere the input is 2.0**100. +c* +c********** +c +c this version just calls zaiord l times. +c + integer k +c + do k= 1,l + call zaiord(h(1,1,k), mask,lmask, + & hmin(k),hmax(k), iaunit) + enddo +c + return + end subroutine zaiord3 + + subroutine zaiord(h, mask,lmask, hmin,hmax, iaunit) + implicit none +c + logical, intent(in) :: lmask + integer, intent(in) :: iaunit + integer, dimension (1:idm,1:jdm), + & intent(in) :: mask +#if defined(REAL4) + real*4, intent(out) :: hmin,hmax + real*4, dimension (1:idm,1:jdm), + & intent(out) :: h +#else + real, intent(out) :: hmin,hmax + real, dimension (1:idm,1:jdm), + & intent(out) :: h +#endif +c +c********** +c* +c 1) machine specific routine for array reading. +c +c must call zaiopn for this array unit before calling zaiord. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c the array, 'h', must conform to that passed in the associated +c call to zaiopn. +c +c 4) hmin,hmax are returned as the minimum and maximum value in the +c array, ignoring array elements set to 2.0**100. +c if lmask==.true. the range is calculated only where mask.ne.0, +c with all other values unchanged in h on exit. It is then an +c error if mask.ne.0 anywhere the input is 2.0**100. +c* +c********** +c + integer ios, i,j + real*4 wmin,wmax +#if defined(TIMER) +c + call xctmr0(17) +#endif +c +* write(6,*) 'zaiord - iaunit,rec = ',iaunit,iarec(iaunit) +* call flush(6) + if (iarec(iaunit).lt.0) then + write(6,9000) iaunit + call flush(6) + stop + endif +c + iarec(iaunit) = iarec(iaunit) + 1 + call zaiordd(w,n2drec, iaunit+1000,iarec(iaunit),ios) + if (ios.ne.0) then + write(6,9100) iarec(iaunit),iaunit + write(6,*) 'ios = ',ios + call flush(6) + stop + endif + wmin = spval + wmax = -spval + if (lmask) then +!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP& REDUCTION(MIN:wmin) REDUCTION(MAX:wmax) + do j= 1,jdm + do i= 1,idm + if (mask(i,j).ne.0) then + h(i,j) = w(i+(j-1)*idm) + wmin = min( wmin, w(i+(j-1)*idm) ) + wmax = max( wmax, w(i+(j-1)*idm) ) + endif + enddo + enddo + if (wmax.eq.spval) then + write(6,9200) iarec(iaunit),iaunit + call flush(6) + stop + endif + else +!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP& REDUCTION(MIN:wmin) REDUCTION(MAX:wmax) + do j= 1,jdm + do i= 1,idm + h(i,j) = w(i+(j-1)*idm) + if (w(i+(j-1)*idm).ne.spval) then + wmin = min( wmin, w(i+(j-1)*idm) ) + wmax = max( wmax, w(i+(j-1)*idm) ) + endif + enddo + enddo + endif + hmin = wmin + hmax = wmax +c +#if defined(TIMER) +c + call xctmr1(17) +#endif + return +c + 9000 format(/ /10x,'error in zaiord - array I/O unit ', + & i3,' is not marked as open.'/ /) + 9100 format(/ /10x,'error in zaiord - can''t read record', + & i4,' on array I/O unit ',i3,'.'/ /) + 9200 format(/ /10x,'error in zaiord - record', + & i4,' on array I/O unit ',i3, + & ' has 2.0**100 outside masked region.'/ /) + end subroutine zaiord + + subroutine zaiordd(a,n, iunit,irec,ios) + implicit none +c + integer, intent(in) :: n,iunit,irec + integer, intent(out) :: ios + real*4, intent(out) :: a(n) +c +c********** +c* +c 1) direct access read a single record. +c +c 2) expressed as a subroutine because i/o with +c implied do loops can be slow on some machines. +c* +c********** +c + read(unit=iunit, rec=irec, iostat=ios) a +#if defined(ENDIAN_IO) + call zaio_endian(a,n) +#endif + return + end subroutine zaiordd + + subroutine zaiosk(iaunit) + implicit none +c + integer, intent(in) :: iaunit +c +c********** +c* +c 1) machine specific routine for skipping an array read. +c +c must call zaiopn for this array unit before calling zaiosk. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c the array, 'h', must conform to that passed in the associated +c call to zaiopn. +c* +c********** +#if defined(TIMER) +c + call xctmr0(16) +#endif +c +* write(6,*) 'zaiosk - iaunit,rec = ',iaunit,iarec(iaunit) +* call flush(6) + if (iarec(iaunit).lt.0) then + write(6,9000) iaunit + call flush(6) + stop + endif +c + iarec(iaunit) = iarec(iaunit) + 1 +#if defined(TIMER) +c + call xctmr1(16) +#endif + return +c + 9000 format(/ /10x,'error in zaiosk - array I/O unit ', + & i3,' is not marked as open.'/ /) + end subroutine zaiosk + + subroutine zaiowr3(h, l, mask,lmask, hmin,hmax, iaunit, lreal4) + implicit none +c + logical, intent(in) :: lmask,lreal4 + integer, intent(in) :: l,iaunit + integer, dimension (1:idm,1:jdm), + & intent(in) :: mask +#if defined(REAL4) + real*4, intent(out) :: hmin(l),hmax(l) + real*4, dimension (1:idm,1:jdm,l), + & intent(inout) :: h +#else + real, intent(out) :: hmin(l),hmax(l) + real, dimension (1:idm,1:jdm,l), + & intent(inout) :: h +#endif +c +c********** +c* +c 1) machine specific routine for 3-d array writing. +c +c must call zaiopn for this array unit before calling zaiord. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c the array, 'h', must conform to that passed in the associated +c call to zaiopn. +c +c 4) hmin,hmax are returned as the minimum and maximum value in the array. +c if lmask==.true. the range is only where mask.ne.0, with all other +c values output as 2.0**100. +c +c 5) If lreal4==.true. then h is overwritten on exit with real*4 version +c of the same array. This is typically used for reproducability on +c restart. +c* +c********** +c +c this version just calls zaiowr l times. +c + integer k +c + do k= 1,l + call zaiowr(h(1,1,k), mask,lmask, + & hmin(k),hmax(k), iaunit, lreal4) + enddo + return + end subroutine zaiowr3 + + subroutine zaiowr(h, mask,lmask, hmin,hmax, iaunit, lreal4) + implicit none +c + logical, intent(in) :: lmask,lreal4 + integer, intent(in) :: iaunit + integer, dimension (1:idm,1:jdm), + & intent(in) :: mask +#if defined(REAL4) + real*4, intent(out) :: hmin,hmax + real*4, dimension (1:idm,1:jdm), + & intent(inout) :: h +#else + real, intent(out) :: hmin,hmax + real, dimension (1:idm,1:jdm), + & intent(inout) :: h +#endif +c +c********** +c* +c 1) machine specific routine for array writing. +c +c must call zaiopn for this array unit before calling zaiord. +c +c 2) array i/o is fortran real*4 direct access i/o to unit iaunit+1000. +c +c 3) iaunit+1000 is the i/o unit used for arrays. array i/o might not +c use fortran i/o units, but, for compatability, assume that +c iaunit+1000 refers to a fortran i/o unit anyway. +c the array, 'h', must conform to that passed in the associated +c call to zaiopn. +c +c 4) hmin,hmax are returned as the minimum and maximum value in the array. +c if lmask==.true. the range is only where mask.ne.0, with all other +c values output as 2.0**100. +c +c 5) If lreal4==.true. then h is overwritten on exit with real*4 version +c of the same array. This is typically used for reproducability on +c restart. +c* +c********** +c + integer ios, i,j + real*4 wmin,wmax +#if defined(TIMER) +c + call xctmr0(18) +#endif +c + if (iarec(iaunit).lt.0) then + write(6,9000) iaunit + call flush(6) + stop + endif +c + wmin = spval + wmax = -spval + if (lreal4) then + if (lmask) then +!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP& REDUCTION(MIN:wmin) REDUCTION(MAX:wmax) + do j= 1,jdm + do i= 1,idm + if (mask(i,j).ne.0) then + w(i+(j-1)*idm) = h(i,j) + wmin = min( wmin, w(i+(j-1)*idm) ) + wmax = max( wmax, w(i+(j-1)*idm) ) + else + w(i+(j-1)*idm) = spval + endif +#if defined(REAL4) +! --- h(i,j) = w(i+(j-1)*idm) ! h is already real*4 +#else + h(i,j) = w(i+(j-1)*idm) ! h is not real*4, so update it +#endif + enddo + enddo + else +!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP& REDUCTION(MIN:wmin) REDUCTION(MAX:wmax) + do j= 1,jdm + do i= 1,idm + w(i+(j-1)*idm) = h(i,j) + if (w(i+(j-1)*idm).ne.spval) then + wmin = min( wmin, w(i+(j-1)*idm) ) + wmax = max( wmax, w(i+(j-1)*idm) ) + endif +#if defined(REAL4) +! --- h(i,j) = w(i+(j-1)*idm) ! h is already real*4 +#else + h(i,j) = w(i+(j-1)*idm) ! h is not real*4, so update it +#endif + enddo + enddo + endif + else + if (lmask) then +!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP& REDUCTION(MIN:wmin) REDUCTION(MAX:wmax) + do j= 1,jdm + do i= 1,idm + if (mask(i,j).ne.0) then + w(i+(j-1)*idm) = h(i,j) + wmin = min( wmin, w(i+(j-1)*idm) ) + wmax = max( wmax, w(i+(j-1)*idm) ) + else + w(i+(j-1)*idm) = spval + endif + enddo + enddo + else +!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP& REDUCTION(MIN:wmin) REDUCTION(MAX:wmax) + do j= 1,jdm + do i= 1,idm + w(i+(j-1)*idm) = h(i,j) + if (w(i+(j-1)*idm).ne.spval) then + wmin = min( wmin, w(i+(j-1)*idm) ) + wmax = max( wmax, w(i+(j-1)*idm) ) + endif + enddo + enddo + endif + endif + do i= idm*jdm+1,n2drec + w(i) = spval + enddo + hmin = wmin + hmax = wmax + iarec(iaunit) = iarec(iaunit) + 1 + call zaiowrd(w,n2drec, iaunit+1000,iarec(iaunit),ios) + if (ios.ne.0) then + write(6,9100) iarec(iaunit),iaunit + call flush(6) + stop + endif +#if defined(TIMER) +c + call xctmr1(18) +#endif + return +c + 9000 format(/ /10x,'error in zaiowr - array I/O unit ', + & i3,' is not marked as open.'/ /) + 9100 format(/ /10x,'error in zaiowr - can''t write record', + & i4,' on array I/O unit ',i3,'.'/ /) + end subroutine zaiowr + + subroutine zaiowrd(a,n, iunit,irec,ios) + implicit none +c + integer, intent(in) :: n,iunit,irec + integer, intent(out) :: ios + real*4, intent(in) :: a(n) +c +c********** +c* +c 1) direct access write a single record. +c +c 2) expressed as a subroutine because i/o with +c implied do loops can be slow on some machines. +c* +c********** +c +#if defined(ENDIAN_IO) + call zaio_endian(a,n) ! overwrites a +#endif + write(unit=iunit, rec=irec, iostat=ios) a + return + end subroutine zaiowrd + + end module mod_za + +#if defined(ENDIAN_IO) + subroutine zaio_endian(a,n) + implicit none +c + integer, intent(in) :: n + integer(kind=4), intent(inout) :: a(n) ! 4-bytes +c +c********** +c* +c 1) swap the endian-ness of the array. +c +c 2) assumes integer(kind=1) and integer(kind=4) ocupy one and four +c bytes respectively. +c* +c********** +c + integer k + integer(kind=4) ii4, io4 ! 4-bytes + integer(kind=1) ii1(4),io1(4) ! 1-byte + equivalence (ii4,ii1(1)), (io4,io1(1)) ! non-standard f90 +c + do k= 1,n + ii4 = a(k) + io1(1) = ii1(4) + io1(2) = ii1(3) + io1(3) = ii1(2) + io1(4) = ii1(1) + a(k) = io4 + enddo + return + end subroutine zaio_endian +#endif /* ENDIAN_IO */ diff --git a/sorc/hycom_tools_for_mom6.fd/ncodaz_inc2mom6nc_glb_lyr.f b/sorc/hycom_tools_for_mom6.fd/ncodaz_inc2mom6nc_glb_lyr.f new file mode 100644 index 00000000..470f3cdd --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/ncodaz_inc2mom6nc_glb_lyr.f @@ -0,0 +1,770 @@ + program ncodaz_inc2mom6nc_glb + use mod_mom6 ! HYCOM mom6 array interface + use mod_za ! HYCOM array I/O interface + use mod_ncio ! A. Srinivasan tsis netcdf module + use netcdf ! Netcdf module + +c +c modified from a COAPS utility by Alex Bozec + implicit none +c +c --- convert NCODA layer increments to MOM6 Netcdf for inc. update +c --- requires MOM6 restart template and +c --- the HYCOM regional.grid for the MOM6 domain. +c --- uses the ncoda increments in the background layers, and the +c --- ncoda binary containing the background state layer thicknesses +c --- which add to depth +c + character*256 flnm_it,flnm_is,flnm_iu,flnm_iv + character*256 flnm_ots,flnm_ouv,filename + character*256 flnm_t,flnm_s,flnm_u,flnm_v,flnm_p,flnm_h + logical larctic,lsymetr + integer i,ia,j,k,ntq,mtq,mro,ib,jb + integer kmax + integer yrflag + real xmin,xmax,misval + real zl + integer itest, jtest + double precision time3(3) +c --- NCODA increments + integer kncoda + real, allocatable :: incoda(:,:) + real, allocatable :: sum(:,:) + + real, allocatable :: tncoda(:,:,:),sncoda(:,:,:) + & ,pncoda(:,:,:) + real, allocatable :: uncoda(:,:,:),vncoda(:,:,:) + real, allocatable :: arrayncoda(:,:,:) + +c --- MOM6 restart fields + integer fid + character(:),allocatable :: vname + integer :: vtype + integer :: vdim1d(1),vdim3d(3),vdim4d(4) +c --- 4D fields + real, allocatable, dimension (:,:,:,:) :: + & u_inc,v_inc,temp_inc,saln_inc + real, allocatable, dimension (:,:,:,:) :: zh +c --- 1D fields + real, allocatable, dimension (:) :: + & lath,lonh,latq,lonq,Layer,intd,time + & ,zz + + + lp = 6 + yrflag = 3 + +c +c --- 'flnm_it' = name of mom6 restart file containing pot. temp. (input) +c --- 'flnm_is' = name of mom6 restart file containing salinity (input) +c --- 'flnm_iu' = name of mom6 restart file containing U (input) +c --- 'flnm_iv' = name of mom6 restart file containing V (input) +c +c --- 'flnm_ots' = name of MOM6 increment file for T,S and h (output) +c --- 'flnm_ouv' = name of MOM6 increment file for U, V and h (output) +c + read (*,'(a)') flnm_it + read (*,'(a)') flnm_is + read (*,'(a)') flnm_iu + read (*,'(a)') flnm_iv + write (lp,'(2a)') ' input MOM6 T file: ',trim(flnm_it) + write (lp,'(2a)') ' input MOM6 S file: ',trim(flnm_is) + write (lp,'(2a)') ' input MOM6 U file: ',trim(flnm_iu) + write (lp,'(2a)') ' input MOM6 V file: ',trim(flnm_iv) + +c read (*,'(a)') flnm_h +c write (lp,'(2a)') 'input ncoda-written background state h file:', +c & trim(flnm_h) +c call flush(lp) + read (*,'(a)') flnm_ots + write (lp,'(2a)') 'output MOM6 TSh file: ',trim(flnm_ots) + read (*,'(a)') flnm_ouv + write (lp,'(2a)') 'output MOM6 UV file: ',trim(flnm_ouv) + if (trim(flnm_ots) .eq. trim(flnm_ouv)) then + write (lp,'(2a)') 'output MOM6 TSh and UV should be different' + call flush(lp) + stop + endif + call flush(lp) +c --- 'ii ' = zonal size of the NCODA increments +c --- 'jj ' = meridional size of the MOM6 ncoda increments +c --- 'kk ' = number of layers in restart +!c --- 'lsymet' = T model is symmetric, F model is non-symmetric +!c --- 'larcti' = T model has an Arctic patch + call blkini(ii,'ii ') + call blkini(jj,'jj ') + call blkini(kk,'kk ') + !call blkinl(lsymetr,'lsymet') + !call blkinl(larctic,'larcti') + + + +c +c --- allocate a few state variables +c + + + + call mom6_alloc +c --- mom6 dimensions +c + call rd_dimen(ntq,mto,kk,mro, flnm_iu,'u') + call rd_dimen(nto,mtq,kk,mro, flnm_iv,'v') + misval = -1.e20 !no missing values in input fields +c + lsymetr = mtq .eq. jj+1 .and. ntq .eq. ii+1 + !larctic = mto .eq. jdm-1 .and. nto .eq. idm +c + write(lp,*) + write(lp,*) 'nto,mto = ',nto,mto + write(lp,*) 'ntq,mtq = ',ntq,mtq + write(lp,*) 'ii,jj,kk = ',ii, jj, kk + write(lp,*) 'lsymetr = ',lsymetr + !write(lp,*) 'larctic = ',larctic + write(lp,*) + write(lp,*) 'misval = ',misval + write(lp,*) + call zhflsh(lp) +c +c +c --- some array allocations are in-line below +c + write(lp,*) 'before Allocate' + + allocate(lath(mto)) + allocate(lonh(nto)) + allocate(latq(mtq)) + allocate(lonq(ntq)) + allocate(Layer(kk)) + + allocate(time(1)) + write(lp,*) 'After allocate' + +c +c --- read the mom6 file. +c +c irec = 1 +c --- read invariant variables + write(lp,*) 'File:',flnm_it + !write(lp,*) 'fid:',fid + +c --- get coordinates t-points + call nciopn(flnm_it,fid) + call nciorv(flnm_it,fid,"lath",lath(:)) + call nciorv(flnm_it,fid,"lonh",lonh(:)) + call nciorv(flnm_it,fid,"Layer",layer(:)) + call nciorv(flnm_it,fid,"Time",time(:)) + call nciocl(flnm_it,fid) +c --- get coordinates u-points + call nciopn(flnm_iu,fid) + call nciorv(flnm_iu,fid,"lonq",lonq(:)) + call nciocl(flnm_iu,fid) +c --- get coordinates v-points + call nciopn(flnm_iv,fid) + call nciorv(flnm_iv,fid,"latq",latq(:)) + call nciocl(flnm_iv,fid) + + +c --- NCODA inc +c --- 'kncoda' = number of ncoda levels +c --- 'flnm_t' = name of ncoda temperature increment file, or "NONE" to exit +c --- 'flnm_s' = name of ncoda salinity increment file +c --- 'flnm_u' = name of ncoda u-velocity increment file, or "NONE" +c --- 'flnm_v' = name of ncoda v-velocity increment file, or "NONE" +c --- 'flnm_p' = name of ncoda density displacement file, or "NONE" +c --- 'flnm_h' = name of the background field layer thick. as ncoda binary + +c itest,jtest show be read with blkini2, now hard coded + itest=1324 + jtest=1898 + call blkini(kncoda,'kncoda') + if (kncoda.ne.kk) then + write(6,*)'error, kncoda has to be kk for lyr option' + stop + endif +c --- allocate ncoda fields + allocate(temp_inc(nto,mto,kncoda,1)) + allocate(saln_inc(nto,mto,kncoda,1)) + allocate(u_inc(ntq,mto,kncoda,1)) + allocate(v_inc(nto,mtq,kncoda,1)) + allocate(zh(nto,mto,kncoda,1)) + +c --- read the nominal layer zz, for example the values in coordinate 'Layer' +c --- in the mom6 archives + allocate (zz(kncoda)) + do k=1,kncoda + call blkinr(zl,'zl ','(a6," =",f10.4)') + zz(k)=zl + enddo + +c --- zh is the background field thickness, which is already normalized to add to depth +c ----read from an ncoda restart file + read (*,'(a)') flnm_t + write (lp,'(2a)') 'Tncoda file: ',trim(flnm_t) + if (flnm_t.eq."NONE") then + write(lp,*) + write(lp,*) '***** EXIT ncoda_inc2mom6ncr *****' + write(lp,*) + call flush(lp) + stop !ncoda_cycle + endif + read (*,'(a)') flnm_s + write (lp,'(2a)') 'Sncoda file: ',trim(flnm_s) + read (*,'(a)') flnm_u + write (lp,'(2a)') 'Uncoda file: ',trim(flnm_u) + read (*,'(a)') flnm_v + write (lp,'(2a)') 'Vncoda file: ',trim(flnm_v) + read (*,'(a)') flnm_p + write (lp,'(2a)') 'Pncoda file: ',trim(flnm_p) + read (*,'(a)') flnm_h + write (lp,'(2a)') 'h background file: ',trim(flnm_p) + +C + +c --- read increments +c --- deallocate these arrays at end of ncoda_cycle loop + allocate( tncoda(ii,jj,kncoda), + & sncoda(ii,jj,kncoda), + & pncoda(ii,jj,kncoda), + & arrayncoda(ii,jj,kncoda)) +c + print*,'Alex idm,jdm,jdmn,kncoda:',ii,jj,kncoda + + write(lp,*) 'open ',trim(flnm_t) + call flush(lp) +c read all EMC ncoda binaries as stream + open(9,file=flnm_t,form='unformatted',access='stream', + & status='old') + + read(unit=9) arrayncoda + do k=1,kncoda + tncoda(:,:,k)=arrayncoda(:,1:jj,k) + enddo +c*** + do k=1,kncoda + write(lp,*) 'Temp incr = ',minval(tncoda(:,:,k)), + & maxval(tncoda(:,:,k)) + if (min(itest,jtest).gt.0) + & write(lp,*) itest,jtest,'Temp incr= ', tncoda(itest,jtest,k) + enddo + + close(unit=9) + write(lp,*) 'close ',trim(flnm_t) + call flush(lp) +c + write(lp,*) 'open ',trim(flnm_s) + call flush(lp) + open(9,file=flnm_s,form='unformatted',access='stream', + & status='old') + read(unit=9) arrayncoda + do k=1,kncoda + sncoda(:,:,k)=arrayncoda(:,1:jj,k) + enddo +c*** + do k=1,kncoda + write(lp,*) 'Salt incr = ',minval(sncoda(:,:,k)), + & maxval(sncoda(:,:,k)) + if (min(itest,jtest).gt.0) + & write(lp,*) 'Salt incr= ', sncoda(itest,jtest,k) + enddo + + close(unit=9) + write(lp,*) 'close ',trim(flnm_s) + call flush(lp) +c + if (flnm_u.ne."NONE") then +c --- deallocate these arrays at end of ncoda_cycle loop + allocate( uncoda(ii,jj,kncoda), + & vncoda(ii,jj,kncoda) ) +c + write(lp,*) 'open ',trim(flnm_u) + call flush(lp) + + + open(9,file=flnm_u,form='unformatted',access='stream', + & status='old') + read(unit=9) arrayncoda + + do k=1,kncoda + do j=1,jj + do i=1,ii +c if (i.eq.1) then +c ib=ii !periodic +c else +c ib=i-1 +c endif +c uncoda(i,j,k)=0.5*(arrayncoda(i,j,k)+arrayncoda(ib,j,k)) +c uncoda is re-staggered in ncoda + uncoda(i,j,k)=arrayncoda(i,j,k) + enddo !i + enddo !j + enddo !k +c*** + do k=1,kncoda + write(lp,*) 'uvel incr = ',minval(uncoda(:,:,k)), + & maxval(uncoda(:,:,k)) + enddo + + close(unit=9) + write(lp,*) 'close ',trim(flnm_u) + call flush(lp) +c + write(lp,*) 'open ',trim(flnm_v) + call flush(lp) + open(9,file=flnm_v,form='unformatted',access='stream', + & status='old') + read(unit=9) arrayncoda + + do k=1,kncoda + do j=1,jj +c jb=max(j-1,1) + do i=1,ii +c vncoda(i,j,k)=0.5*(arrayncoda(i,j,k)+arrayncoda(i,jb,k)) +c vncoda is re-staggered in ncoda + vncoda(i,j,k)=arrayncoda(i,j,k) + enddo !i + enddo !j + enddo !k +c*** + do k=1,kncoda + write(lp,*) 'vvel incr = ',minval(vncoda(:,:,k)), + & maxval(vncoda(:,:,k)) + enddo + + close(unit=9) + write(lp,*) 'close ',trim(flnm_v) + call flush(lp) + endif !u&vncoda +c + if (flnm_p.ne."NONE") then + write(lp,*) 'open ',trim(flnm_p) + call flush(lp) + open(9,file=flnm_p,form='unformatted',access='stream', + & status='old') + read(unit=9) arrayncoda + do k=1,kncoda + pncoda(:,:,k)=arrayncoda(:,1:jj,k) + enddo +c*** + do k=1,kncoda + write(lp,*) 'intd incr = ',minval(pncoda(:,:,k)), + & maxval(pncoda(:,:,k)) + enddo + close(unit=9) + write(lp,*) 'close ',trim(flnm_p) + call flush(lp) + endif !pncoda + + +c --- get MOM6 z-layers from an ncoda restart + open(9, file=flnm_h, form='unformatted', access='stream', + & status='old') +c --- ZG sum h has to add to depth, already does + read(9)arrayncoda + do k=1,kncoda + do j= 1,jj + do i= 1,ii + zh(i,j,k,1)=arrayncoda(i,j,k) + enddo + enddo + write(lp,*) 'zh backgr = ',minval(zh(:,:,k,1)), + & maxval(zh(:,:,k,1)) + enddo + +c --- store increments in proper dimension arrays + temp_inc(:,:,:,1) = 0.0 + saln_inc(:,:,:,1) = 0.0 + do j= 1,mto + do i= 1,nto + temp_inc(i,j,:,1)=tncoda(i,j,:) + saln_inc(i,j,:,1)=sncoda(i,j,:) + enddo + enddo + + u_inc(:,:,:,1) = 0.0 + do j= 1,mto + do i= 1,nto + u_inc(i,j,:,1)=uncoda(i,j,:) + enddo + enddo +c --- fill last column + if (lsymetr) then + do j= 1,mto + i=ntq + u_inc(i,j,:,1)=u_inc(1,j,:,1) + enddo + endif + + v_inc(:,:,:,1) = 0.0 + do j= 1,mto + do i= 1,nto + v_inc(i,j,:,1)=vncoda(i,j,:) + enddo + enddo +c --- fill last row + if (lsymetr) then +!!Alex not do anything for now since no velocity correction in Arctic +!NOTE in the future see MOM6 handles Arctic patch and reproduce for v +! do i= 1,nto +! j=mtq +! v_inc(i,j,:,1)=v_inc(i,mtq-1,:,1) +! enddo + endif +c --- write MOM6 TS,h increment file + !create file + filename=flnm_ots + + CALL nciocf(filename,fid) + + !create dimension and coordinate variable 1 + vname="lonh" + vtype=NF90_DOUBLE + vdim1d(1)=1 + CALL nciodd(filename,fid,vname,nto) + CALL nciocv(filename,fid,vname,vtype,vdim1d) + CALL nciowa(filename,fid,vname,'long_name','Longitude') + CALL nciowa(filename,fid,vname,'cartesian_axis','X') + CALL nciowa(filename,fid,vname,'units','degrees_east') + + + !create dimension and coordinate variable 2 + vname="lath" + vtype=NF90_DOUBLE + vdim1d(1)=2 + CALL nciodd(filename,fid,vname,mto) + CALL nciocv(filename,fid,vname,vtype,vdim1d) + CALL nciowa(filename,fid,vname,'long_name','Latitude') + CALL nciowa(filename,fid,vname,'cartesian_axis','Y') + CALL nciowa(filename,fid,vname,'units','degrees_north') + + !create dimension and coordinate variable 5 +c vname="Level" + vname="Layer" + vtype=NF90_DOUBLE + vdim1d(1)=3 + CALL nciodd(filename,fid,vname,kncoda) +c CALL nciodd(filename,fid,vname,kncoda+1) + CALL nciocv(filename,fid,vname,vtype,vdim1d) +c CALL nciowa(filename,fid,vname,'long_name','z-levels') + CALL nciowa(filename,fid,vname,'long_name','Layer z-rho') + CALL nciowa(filename,fid,vname,'cartesian_axis','Z') + CALL nciowa(filename,fid,vname,'units','meter') + CALL nciowa(filename,fid,vname,'positive','up') + print*,'after Dim' + +! vname="Time" +! vtype=NF90_DOUBLE +! vdim1d(1)=6 +! CALL nciodd(filename,fid,vname,NF90_UNLIMITED) +! CALL nciocv(filename,fid,vname,vtype,vdim1d) +! CALL nciowa(filename,fid,vname,'long_name','Time') +! CALL nciowa(filename,fid,vname,'cartesian_axis','T') +! CALL nciowa(filename,fid,vname,'units','days') + + + !create data variable + vname="pt_inc" + vtype=NF90_DOUBLE + vdim3d=(/1,2,3/) + call nciocv(filename,fid,vname,vtype,vdim3d) + call nciowa(filename,fid,vname,'long_name', + & 'NCODA Potential Temperature increments') + call nciowa(filename,fid,vname,'units','degC') + + !create data variable + vname="s_inc" + vtype=NF90_DOUBLE + vdim3d=(/1,2,3/) + call nciocv(filename,fid,vname,vtype,vdim3d) + call nciowa(filename,fid,vname,'long_name', + & 'NCODA Salinity increments') + call nciowa(filename,fid,vname,'units','PPT') + + !create data variable + vname="zh" + vtype=NF90_DOUBLE + vdim3d=(/1,2,3/) + call nciocv(filename,fid,vname,vtype,vdim3d) + call nciowa(filename,fid,vname,'long_name', + & 'MOM6 z-levels thickness') + call nciowa(filename,fid,vname,'units','m') + + !end define mode + call ncioed(filename,fid) + print*, 'after var defnition' + + ! write data into variables + call nciowv(filename,fid,"lath",lath) + call nciowv(filename,fid,"lonh",lonh) + call nciowv(filename,fid,"Layer",zz(1:kncoda)) + print*, 'after writing Layer' + !call nciowv(filename,fid,"Time",time) + print*, 'after writing dim' + + call nciowv(filename,fid,"pt_inc",temp_inc) + write(lp,*) 'Pot. temperature inc. written' + call nciowv(filename,fid,"s_inc",saln_inc) + write(lp,*) 'Salinity inc. written' + call nciowv(filename,fid,"zh",zh) + write(lp,*) 'backgrnd thk. written' +C + ! close file + call nciocl(filename,fid) + +c --- write MOM6 U,V,h increment file + !create file + filename=flnm_ouv + + CALL nciocf(filename,fid) + + !create dimension and coordinate variable 1 + vname="lonh" + vtype=NF90_DOUBLE + vdim1d(1)=1 + CALL nciodd(filename,fid,vname,nto) + CALL nciocv(filename,fid,vname,vtype,vdim1d) + CALL nciowa(filename,fid,vname,'long_name','Longitude') + CALL nciowa(filename,fid,vname,'cartesian_axis','X') + CALL nciowa(filename,fid,vname,'units','degrees_east') + + + !create dimension and coordinate variable 2 + vname="lath" + vtype=NF90_DOUBLE + vdim1d(1)=2 + CALL nciodd(filename,fid,vname,mto) + CALL nciocv(filename,fid,vname,vtype,vdim1d) + CALL nciowa(filename,fid,vname,'long_name','Latitude') + CALL nciowa(filename,fid,vname,'cartesian_axis','Y') + CALL nciowa(filename,fid,vname,'units','degrees_north') + + !create dimension and coordinate variable 5 + vname="Level" + vtype=NF90_DOUBLE + !vdim1d(1)=5 + vdim1d(1)=3 + CALL nciodd(filename,fid,vname,kncoda) + CALL nciocv(filename,fid,vname,vtype,vdim1d) + CALL nciowa(filename,fid,vname,'long_name','z-levels') + CALL nciowa(filename,fid,vname,'cartesian_axis','Z') + CALL nciowa(filename,fid,vname,'units','meter') + CALL nciowa(filename,fid,vname,'positive','up') + print*,'after Dim' + + !create dimension and coordinate variable 3 + vname="lonq" + vtype=NF90_DOUBLE + !vdim1d(1)=3 + vdim1d(1)=4 + CALL nciodd(filename,fid,vname,ntq) + CALL nciocv(filename,fid,vname,vtype,vdim1d) + CALL nciowa(filename,fid,vname,'long_name','Longitude') + CALL nciowa(filename,fid,vname,'cartesian_axis','X') + CALL nciowa(filename,fid,vname,'units','degrees_east') + + !create dimension and coordinate variable 4 + vname="latq" + vtype=NF90_DOUBLE + !vdim1d(1)=4 + vdim1d(1)=5 + CALL nciodd(filename,fid,vname,mtq) + CALL nciocv(filename,fid,vname,vtype,vdim1d) + CALL nciowa(filename,fid,vname,'long_name','Latitude') + CALL nciowa(filename,fid,vname,'cartesian_axis','Y') + CALL nciowa(filename,fid,vname,'units','degrees_north') + +! vname="Time" +! vtype=NF90_DOUBLE +! vdim1d(1)=6 +! CALL nciodd(filename,fid,vname,NF90_UNLIMITED) +! CALL nciocv(filename,fid,vname,vtype,vdim1d) +! CALL nciowa(filename,fid,vname,'long_name','Time') +! CALL nciowa(filename,fid,vname,'cartesian_axis','T') +! CALL nciowa(filename,fid,vname,'units','days') + + + !create data variable + vname="u_inc" + vtype=NF90_DOUBLE + vdim3d=(/4,2,3/) + call nciocv(filename,fid,vname,vtype,vdim3d) + call nciowa(filename,fid,vname,'long_name', + & 'NCODA Zonal velocity increments') + call nciowa(filename,fid,vname,'units','m s-1') + + !create data variable + vname="v_inc" + vtype=NF90_DOUBLE + vdim3d=(/1,5,3/) + call nciocv(filename,fid,vname,vtype,vdim3d) + call nciowa(filename,fid,vname,'long_name', + & 'NCODA Meridional velocity increments') + call nciowa(filename,fid,vname,'units','m s-1') + + !end define mode + call ncioed(filename,fid) + print*, 'after var defnition' + + ! write data into variables + call nciowv(filename,fid,"lath",lath) + call nciowv(filename,fid,"lonh",lonh) + call nciowv(filename,fid,"latq",latq) + call nciowv(filename,fid,"lonq",lonq) + call nciowv(filename,fid,"Level",zz(1:kncoda)) + !call nciowv(filename,fid,"Time",time) + print*, 'after writing dim' + + call nciowv(filename,fid,"u_inc",u_inc) + write(lp,*) 'U-vel inc. written' + call nciowv(filename,fid,"v_inc",v_inc) + write(lp,*) 'V-vel inc. written' +C + ! close file U,V,h + call nciocl(filename,fid) + + end program ncodaz_inc2mom6nc_glb + + subroutine m2h_p(f_nc,nto,mto,kk,misval, + & field,ii,jj, lsymetr,larctic) + implicit none +c + logical lsymetr,larctic + integer nto,mto,kk,ii,jj + real f_nc(nto,mto,kk),misval,field(ii,jj,kk) +c +c --- spval = hycom data void marker, 2^100 or about 1.2676506e30 +******real, parameter :: spval=2.0**100 + real, parameter :: spval=0.0 !no spval, use 0.0 +c +c --- convert p-grid mom6 array to hycom. +c + integer i,ia,j,k +c + do k= 1,kk + do j= 1,mto + do i= 1,nto + if (f_nc(i,j,k).ne.misval) then + field(i,j,k) = f_nc(i,j,k) + else + field(i,j,k) = spval + endif + enddo !i + enddo !j + if (lsymetr) then + do i= 1,nto + field(i,jj,k) = spval + enddo !i + do j= 1,jj + field(ii,j,k) = spval + enddo !j + elseif (larctic) then !p-grid scalar field, mto=jj-1 + do i= 1,nto + ia = nto-mod(i-1,nto) + field(i,jj,k) = field(ia,jj-1,k) + enddo !i + endif !lsymetr:larctic + enddo !k + return + end + + subroutine m2h_u(f_nc,ntq,mto,kk,misval, + & field,ii,jj,lsymetr,larctic) + implicit none +c + logical lsymetr,larctic + integer ntq,mto,kk,ii,jj + real f_nc(ntq,mto,kk),misval,field(ii,jj,kk) +c +c --- convert u-grid mom6 array to hycom u-grid. +c --- mom6 standard has "q" at i+0.5,j+0.5 w.r.t. p.ij +c --- mom6 symetric has "q" at i-0.5,j-0.5 w.r.t. p.ij +c --- hycom has "q" at i-0.5,j-0.5 w.r.t. p.ij +c + integer i,ia,j,k +c + if (lsymetr) then + do k= 1,kk + do j= 1,mto + do i= 1,ntq + if (f_nc(i,j,k).ne.misval) then + field(i,j,k) = f_nc(i,j,k) + else + field(i,j,k) = 0.0 + endif + enddo !i + enddo !j + do i= 1,ii !must be land, i.e. zero + field(i,jj,k) = 0.0 + enddo !i + enddo !k + else + do k= 1,kk + do j= 1,mto + do i= 1,ntq + ia = mod(i,ntq)+1 + if (f_nc(i,j,k).ne.misval) then + field(ia,j,k) = f_nc(i,j,k) + else + field(ia,j,k) = 0.0 + endif + enddo !i + enddo !j + if (larctic) then !u-grid vector field, mto=jj-1 + do i= 1,ntq + ia = mod(ntq-(i-1),ntq)+1 + field(i,jj,k) = -field(ia,jj-1,k) + enddo !i + endif + enddo !k + endif !lsymetr:else + return + end + + subroutine m2h_v(f_nc,nto,mtq,kk,misval, + & field,ii,jj,lsymetr,larctic) + implicit none +c + logical lsymetr,larctic + integer nto,mtq,kk,ii,jj + real f_nc(nto,mtq,kk),misval,field(ii,jj,kk) +c +c --- convert v-grid mom6 array to hycom. +c --- mom6 standard has "q" at i+0.5,j+0.5 w.r.t. p.ij +c --- mom6 symetric has "q" at i-0.5,j-0.5 w.r.t. p.ij +c --- hycom has "q" at i-0.5,j-0.5 w.r.t. p.ij +c + integer i,ia,j,k +c + if (lsymetr) then + do k= 1,kk + do j= 1,mtq + do i= 1,nto + if (f_nc(i,j,k).ne.misval) then + field(i,j,k) = f_nc(i,j,k) + else + field(i,j,k) = 0.0 + endif + enddo !i + enddo !j + do j= 1,jj !must be land + field(ii,j,k) = 0.0 + enddo !j + enddo !k + else + do k= 1,kk + do j= 1,min(mtq,jj-1) !mtq if larctic + do i= 1,nto + if (f_nc(i,j, k).ne.misval) then + field(i,j+1,k) = f_nc(i,j,k) + else + field(i,j+1,k) = 0.0 + endif + enddo !i + enddo !j + do i= 1,nto !must be land + field(i,1,k) = 0.0 + enddo !i + enddo !k + endif !lsymetr:else + return + end diff --git a/sorc/hycom_tools_for_mom6.fd/putdat.f b/sorc/hycom_tools_for_mom6.fd/putdat.f new file mode 100644 index 00000000..49130589 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/putdat.f @@ -0,0 +1,394 @@ + subroutine putdat(flnm, artype,time, + & lsteric,icegln,trcout, + & iexpt,iversn,yrflag,kkout, thbase) + use mod_mom6 ! HYCOM mom6 array interface + use mod_za ! HYCOM array I/O interface +c + character flnm*(*) + double precision time(3) + real thbase + integer artype,iexpt,iversn,yrflag,kkout + logical lsteric,icegln,trcout +c +c --- write model fields. +c --- HYCOM 2.0 array I/O archive file. +c + character*8 ctype +c + real coord,xmin,xmax + integer i,j,k,ktr,jversn,l,nop + data nop/24/ +c +c +c --- output is in "*.[ab]" +c + l = len_trim(flnm) + if (flnm(l-1:l).eq.'.a' .or. flnm(l-1:l).eq.'.b') then + open (unit=nop,file=flnm(1:l-2)//'.b',form='formatted', + & status='new',action='write') + call zaiopf(flnm(1:l-2)//'.a','new', nop) + else + open (unit=nop,file=flnm(1:l)//'.b',form='formatted', + & status='new',action='write') + call zaiopf(flnm(1:l)//'.a','new', nop) + endif +c +c --- header. +c + jversn = max(iversn,20) + if (artype.eq.1) then + write(nop,116) ctitle,jversn,iexpt,yrflag,idm,jdm + write( lp, *) + write( lp,116) ctitle,jversn,iexpt,yrflag,idm,jdm + 116 format (a80/a80/a80/a80/ + & i5,4x,'''iversn'' = hycom version number x10'/ + & i5,4x,'''iexpt '' = experiment number x10'/ + & i5,4x,'''yrflag'' = days in year flag'/ + & i5,4x,'''idm '' = longitudinal array size'/ + & i5,4x,'''jdm '' = latitudinal array size'/ + & 'field time step model day', + & ' k dens min max') + elseif (artype.eq.2) then + write(nop,118) ctitle,jversn,iexpt,yrflag,idm,jdm + write( lp, *) + write( lp,118) ctitle,jversn,iexpt,yrflag,idm,jdm + 118 format (a80/a80/a80/a80/ + & i5,4x,'''iversn'' = hycom version number x10'/ + & i5,4x,'''iexpt '' = experiment number x10'/ + & i5,4x,'''yrflag'' = days in year flag'/ + & i5,4x,'''idm '' = longitudinal array size'/ + & i5,4x,'''jdm '' = latitudinal array size'/ + & 'field no. recs mean day', + & ' k dens min max') + else + write( lp,"(/ a /)") + & 'error in putdat - only artpe==1 and artype==2 allowed' + stop + endif + +c +c --- surface fields. +c + coord=0.0 +c + call zaiowr(montg,ip,.true., + & xmin,xmax, nop, .false.) + if (sigver.gt.0) then +c --- identify the equation of state on the first record + write (nop,117) 'montg1 ',nstep,time(1),sigver,thbase,xmin,xmax + call flush(nop) + write ( lp,117) 'montg1 ',nstep,time(1),sigver,thbase,xmin,xmax + call flush( lp) + else + write (nop,117) 'montg1 ',nstep,time(1),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'montg1 ',nstep,time(1),0,coord,xmin,xmax + call flush( lp) + endif !sigver + call zaiowr(srfht,ip,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'srfhgt ',nstep,time(2),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'srfhgt ',nstep,time(2),0,coord,xmin,xmax + call flush( lp) +c + if(lsteric) then + call zaiowr(steric,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'steric ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'steric ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + endif +c + call zaiowr(surflx,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'surflx ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'surflx ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(salflx,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'salflx ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'salflx ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) +c + call zaiowr(dpbl,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'bl_dpth ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'bl_dpth ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(dpmixl,ip,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'mix_dpth',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'mix_dpth',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + if (sigver.eq.0) then + call zaiowr(tmix,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'tmix ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'tmix ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(smix,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'smix ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'smix ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(thmix,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'thmix ',nstep,time(3),0,thbase,xmin,xmax + call flush(nop) + write ( lp,117) 'thmix ',nstep,time(3),0,thbase,xmin,xmax + call flush( lp) + call zaiowr(umix,iu,.true., xmin,xmax, nop, .false.) + write (nop,117) 'umix ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'umix ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(vmix,iv,.true., xmin,xmax, nop, .false.) + write (nop,117) 'vmix ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'vmix ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + endif !sigver==0 + if(artype.gt.1) then + call zaiowr(kemix, ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'kemix ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'kemix ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + endif + if(icegln) then + call zaiowr(covice,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'covice ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'covice ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(thkice,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'thkice ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'thkice ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(temice,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'temice ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'temice ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + endif +c +c --- depth averaged fields (no mask, in case this is a subregion archive). +c + call zaiowr(ubaro,iu,.false., + & xmin,xmax, nop, .false.) + write (nop,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(vbaro,iv,.false., + & xmin,xmax, nop, .false.) + write (nop,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + if(artype.gt.1) then + call zaiowr(kebaro,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'kebtrop ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'kebtrop ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + endif +c +c --- layer loop. +c + do 75 k=1,kkout + coord=theta(k) + call zaiowr(u(1,1,k),iu,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'u-vel. ',nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'u-vel. ',nstep,time(3),k,coord,xmin,xmax + call flush( lp) + call zaiowr(v(1,1,k),iv,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'v-vel. ',nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'v-vel. ',nstep,time(3),k,coord,xmin,xmax + call flush( lp) + if(artype.gt.1) then + call zaiowr(ke(1,1,k),ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'k.e. ',nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'k.e. ',nstep,time(3),k,coord,xmin,xmax + call flush( lp) + endif + call zaiowr(dp(1,1,k), ip,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'thknss ',nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'thknss ',nstep,time(3),k,coord,xmin,xmax + call flush( lp) + call zaiowr(temp(1,1,k),ip,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'temp ',nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'temp ',nstep,time(3),k,coord,xmin,xmax + call flush( lp) + call zaiowr(saln(1,1,k),ip,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'salin ',nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'salin ',nstep,time(3),k,coord,xmin,xmax + call flush( lp) + if (sigver.eq.0) then + call zaiowr(th3d(1,1,k),ip,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'density ',nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'density ',nstep,time(3),k,coord,xmin,xmax + call flush( lp) + endif !sigver==0 + if (trcout) then + do ktr= 1,ntracr + call zaiowr(trcr(1,1,k,ktr),ip,.true., + & xmin,xmax, nop, .false.) + if (ktr.le.99) then +c --- older postprocessing code will not recognise this + write(ctype,'(a6,i2.2)') 'tracer',ktr + else + ctype = 'tracer ' + endif + write (nop,117) ctype,nstep,time(3),k,coord,xmin,xmax + call flush(nop) + write ( lp,117) ctype,nstep,time(3),k,coord,xmin,xmax + call flush( lp) + enddo !ktr + endif !trcout + 75 continue +c + 117 format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7) +c + close (unit=nop) + call zaiocl(nop) + return + end + + subroutine putdat_tide(flnm, artype,time,lsteric, + & iexpt,iversn,yrflag, thbase) + use mod_mom6 ! HYCOM mom6 array interface + use mod_za ! HYCOM array I/O interface +c + character flnm*(*) + double precision time(3) + real thbase + integer artype,iexpt,iversn,yrflag + logical lsteric +c +c --- write model fields. +c --- HYCOM 2.0 array I/O archive file. +c + character*8 ctype +c + real coord,xmin,xmax + integer i,j,k,ktr,jversn,l,nop + data nop/24/ +c +c +c --- output is in "*.[ab]" +c + l = len_trim(flnm) + if (flnm(l-1:l).eq.'.a' .or. flnm(l-1:l).eq.'.b') then + open (unit=nop,file=flnm(1:l-2)//'.b',form='formatted', + & status='new',action='write') + call zaiopf(flnm(1:l-2)//'.a','new', nop) + else + open (unit=nop,file=flnm(1:l)//'.b',form='formatted', + & status='new',action='write') + call zaiopf(flnm(1:l)//'.a','new', nop) + endif +c +c --- header. +c + jversn = max(iversn,20) + if (artype.eq.1) then + write(nop,116) ctitle,jversn,iexpt,yrflag,idm,jdm + write( lp, *) + write( lp,116) ctitle,jversn,iexpt,yrflag,idm,jdm + 116 format (a80/a80/a80/a80/ + & i5,4x,'''iversn'' = hycom version number x10'/ + & i5,4x,'''iexpt '' = experiment number x10'/ + & i5,4x,'''yrflag'' = days in year flag'/ + & i5,4x,'''idm '' = longitudinal array size'/ + & i5,4x,'''jdm '' = latitudinal array size'/ + & 'field time step model day', + & ' k dens min max') + elseif (artype.eq.2) then + write(nop,118) ctitle,jversn,iexpt,yrflag,idm,jdm + write( lp, *) + write( lp,118) ctitle,jversn,iexpt,yrflag,idm,jdm + 118 format (a80/a80/a80/a80/ + & i5,4x,'''iversn'' = hycom version number x10'/ + & i5,4x,'''iexpt '' = experiment number x10'/ + & i5,4x,'''yrflag'' = days in year flag'/ + & i5,4x,'''idm '' = longitudinal array size'/ + & i5,4x,'''jdm '' = latitudinal array size'/ + & 'field no. recs mean day', + & ' k dens min max') + else + write( lp,"(/ a /)") + & 'error in putdat - only artpe==1 and artype==2 allowed' + stop + endif + +c +c --- surface fields. +c + coord=0.0 +c + call zaiowr(montg,ip,.true., + & xmin,xmax, nop, .false.) + if (sigver.gt.0) then +c --- identify the equation of state on the first record + write (nop,117) 'montg1 ',nstep,time(1),sigver,thbase,xmin,xmax + call flush(nop) + write ( lp,117) 'montg1 ',nstep,time(1),sigver,thbase,xmin,xmax + call flush( lp) + else + write (nop,117) 'montg1 ',nstep,time(1),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'montg1 ',nstep,time(1),0,coord,xmin,xmax + call flush( lp) + endif !sigver + call zaiowr(srfht,ip,.true., + & xmin,xmax, nop, .false.) + write (nop,117) 'srfhgt ',nstep,time(2),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'srfhgt ',nstep,time(2),0,coord,xmin,xmax + call flush( lp) +c + if(lsteric) then + call zaiowr(steric,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'steric ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'steric ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + endif +c +c --- depth averaged fields (no mask, in case this is a subregion archive). +c + call zaiowr(ubaro,iu,.false., + & xmin,xmax, nop, .false.) + write (nop,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) + call zaiowr(vbaro,iv,.false., + & xmin,xmax, nop, .false.) + write (nop,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax + call flush( lp) +c + 117 format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7) +c + close (unit=nop) + call zaiocl(nop) + return + end diff --git a/sorc/hycom_tools_for_mom6.fd/wtime.F b/sorc/hycom_tools_for_mom6.fd/wtime.F new file mode 100644 index 00000000..5617f78c --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/wtime.F @@ -0,0 +1,77 @@ +#if defined(MPI) + REAL*8 FUNCTION WTIME() + IMPLICIT NONE +C +C USE THE MPI FUNCTION MPI_WTIME TO RETURN WALL TIME. +C + DOUBLE PRECISION MPI_WTIME +C + WTIME = MPI_WTIME() + RETURN +C END OF WTIME. + END +#else + REAL*8 FUNCTION WTIME() + IMPLICIT NONE +C +C USE THE F90 INTRINSIC SYSTEM_CLOCK TO RETURN WALL TIME. +C +C WILL FAIL IF THE COUNT IS EVER NEGATIVE, BUT THE STANDARD +C SAYS THAT IT IS AWAYS NON-NEGATIVE IF A CLOCK EXISTS. +C NOT THREAD-SAFE, UNLESS LCOUNT AND IOVER ARE THREADPRIVATE. +C + REAL*8 ZERO,ONE + PARAMETER (ZERO=0.0, ONE=1.0) +C + INTEGER COUNT, MCOUNT, RATE +C + REAL*8 OFFSEC, OFFSET, PERSEC + INTEGER ICOUNT, IOVER, LCOUNT, NCOUNT + SAVE OFFSEC, OFFSET, PERSEC + SAVE ICOUNT, IOVER, LCOUNT, NCOUNT +C + DATA IOVER, LCOUNT / -1, -1 / +C + CALL SYSTEM_CLOCK(COUNT) +C + IF (COUNT.LT.LCOUNT) THEN +C +C COUNT IS SUPPOSED TO BE NON-DECREASING EXCEPT WHEN IT WRAPS, +C BUT SOME IMPLEMENTATIONS DON''T DO THIS. SO IGNORE ANY +C DECREASE OF LESS THAN ONE PERCENT OF THE RANGE. +C + IF (LCOUNT-COUNT.LT.NCOUNT) THEN + COUNT = LCOUNT + ELSE + IOVER = IOVER + 1 + OFFSET = OFFSET + OFFSEC + ENDIF + ENDIF + LCOUNT = COUNT +C + IF (IOVER.EQ.0) THEN +C +C FIRST CYCLE, FOR ACCURACY WITH 64-BIT COUNTS. +C + WTIME = (COUNT - ICOUNT) * PERSEC + ELSEIF (IOVER.GT.0) THEN +C +C ALL OTHER CYCLES. +C + WTIME = COUNT * PERSEC + OFFSET + ELSE +C +C INITIALIZATION. +C + CALL SYSTEM_CLOCK(ICOUNT, RATE, MCOUNT) + NCOUNT = MCOUNT/100 + PERSEC = ONE/RATE + OFFSEC = MCOUNT * PERSEC + OFFSET = -ICOUNT * PERSEC + IOVER = 0 + WTIME = ZERO + ENDIF + RETURN +C END OF WTIME. + END +#endif /* MPI:else */ diff --git a/sorc/hycom_tools_for_mom6.fd/zh.F b/sorc/hycom_tools_for_mom6.fd/zh.F new file mode 100644 index 00000000..58a6f760 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/zh.F @@ -0,0 +1,67 @@ +#if defined(T3E) +# include "zh_t3e.f" +#elif defined(YMP) +# include "zh_ymp.f" +#elif defined(X1) +# include "zh_x1.f" +#elif defined(HPUX) +# include "zh_hp.f" +#else +# include "zh_sun.f" +#endif +# if defined(AIX) + subroutine flush(iunit) + implicit none + integer iunit +c +c --- wrapper for flush system call under AIX. +c + call flush_(iunit) + return + end +#endif +# if defined(OSX) + subroutine flush(iunit) + implicit none + integer iunit +c +c --- wrapper for flush system call under OSX. +c + call flush_(iunit) + return + end + + subroutine fdate(stm) + implicit none + character*24 stm +c +c --- wrapper for fdate (SunOS) system call under OSX. +c --- An example of SunOS fdate output is: "Mon Aug 1 09:24:21 1994". +c + character*26 str +c + call fdate_(str) + stm = str(1:24) + return + end + + real*4 function etime(time) + implicit none + real*4 time(2) +c +c --- wrapper for etime (SunOS) function system call under OSX. +c + real(4) etime_ + type tb_type + sequence + real(4) usrtime + real(4) systime + end type + type (tb_type) etime_struct +c + etime = etime_(etime_struct) + time(1) = etime_struct%usrtime + time(2) = etime_struct%systime + return + end +#endif diff --git a/sorc/hycom_tools_for_mom6.fd/zh_sun.f b/sorc/hycom_tools_for_mom6.fd/zh_sun.f new file mode 100644 index 00000000..b95c5168 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/zh_sun.f @@ -0,0 +1,308 @@ +C +C----------------------------------------------------------------------- +C +C MACHINE DEPENDENT ROUTINES. +C +C THESE ARE FOR THE SUN. +C +C----------------------------------------------------------------------- +C + SUBROUTINE ZHFLSH(IUNIT) + IMPLICIT NONE +C + INTEGER IUNIT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT FLUSHES THE OUTPUT BUFFERS OF +C LOGICAL UNIT 'IUNIT'. +C +C 2) USE ZAIOFL TO FLUSH ARRAY I/O. +C +C 3) THIS VERSION IS FOR SUN WORKSTATIONS. +C IT USES THE 'FLUSH' FORTRAN SYSTEM ROUTINE. +C +C 4) ALAN J. WALLCRAFT, SEPTEMBER 1989. +C* +C********** +C + CALL FLUSH(IUNIT) + RETURN +C END OF ZHFLSH. + END + SUBROUTINE ZHOPEN(IUNIT,CFORM,CSTAT,IRLEN) + IMPLICIT NONE +C + INTEGER IUNIT,IRLEN + CHARACTER*(*) CFORM,CSTAT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE FOR SIMPLE OPEN STATEMENTS. +C +C SEE ALSO, ZHOPNC. +C +C 2) THIS VERSION IS FOR THE SUN FOR SUN FORTRAN. +C THE FILENAME IS TAKEN FROM THE ENVIRONMENT VARIABLE FOR0xx, +C WHERE xx = IUNIT, WITH DEFAULT fort.xx. +C +C 3) CSTAT CAN BE 'SCRATCH', 'OLD', 'NEW', OR 'UNKNOWN'. +C CFORM CAN BE 'FORMATTED' OR 'UNFORMATTED'. +C IRLEN CAN BE ZERO (FOR SEQUENTIAL ACCESS), OR NON-ZERO (FOR DIRECT +C ACCESS INDICATING RECORD LENGTH IN TERMS OF REAL VARIABLES). +C IF IRLEN IS NEGATIVE, THE OUTPUT WILL BE IN IEEE BINARY, IF THAT +C CAPABILITY EXISTS USING STANDARD FORTRAN I/O. THIS CAPABILITY +C IS PRIMARILY TARGETED TO CRAYS, ON OTHER MACHINES -LEN AND LEN +C ARE LIKELY TO DO THE SAME THING. +C +C ON THE SUN, LEN AND -LEN BOTH GIVE IEEE FILES. +C +C 4) FOR F90 COMPILERS, DELIM='QUOTE' IS INCLUDED IN THE OPEN +C STATEMENT WHERE APPROPRIATE. NOTE THAT THE FOLLOWING CALL +C CALL ZHOPEN(6,'FORMATTED','UNKNOWN',0) +C IS LEGAL, AND WOULD HAVE THE EFFECT OF SETTING DELIM='QUOTE' +C FOR STDOUT. IUNIT=6 IS TYPICALLY TREATED AS A SPECIAL CASE. +C ADDITIONALLY, FOR F90 COMPILERS: +C STATUS='NEW' IMPLIES ACTION='WRITE' +C STATUS='OLD' IMPLIES ACTION='READ' +C STATUS='SCRATCH' IMPLIES ACTION='READWRITE' +C +C 5) ALAN J. WALLCRAFT, DECEMBER 1991 AND AUGUST 1993. +C* +C********** +C + INTEGER IOS,NRECL + CHARACTER CFILE*240,CENV*6 +C +C GET FILENAME. +C + WRITE(CENV,1000) IUNIT + CFILE = ' ' + CALL GETENV(CENV,CFILE) + IF (CFILE.EQ.' ') THEN + WRITE(CFILE,1100) IUNIT + ENDIF +C +C OPEN FILE. +C + IF (IRLEN.EQ.0) THEN +C +C SEQUENTIAL. +C + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + IOSTAT=IOS) + ELSEIF (IUNIT.NE.6) THEN + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + IOSTAT=IOS) + ELSE + IOS = 0 + ENDIF + ELSE +C +C UNFORMATTED DIRECT ACCESS. +C + IF (CFORM.NE.'UNFORMATTED' .AND. + + CFORM.NE.'unformatted' ) THEN + WRITE(6,9100) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF +C + IF (IRLEN.LT.0) THEN +C +C IEEE I/O. +C + NRECL = -4*IRLEN + ELSE + NRECL = 4*IRLEN + ENDIF + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ENDIF + ENDIF + IF (IOS.NE.0) THEN + WRITE(6,9000) IUNIT + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + RETURN +C + 1000 FORMAT('FOR',I3.3) + 1100 FORMAT('fort.',I2.2) + 9000 FORMAT(// 10X,'ERROR IN ZHOPEN - CAN''T OPEN UNIT',I3,'.' //) + 9100 FORMAT(// 10X,'ERROR IN ZHOPEN (UNIT',I3.2,') -' / + + 20X,'ONLY UNFORMATTED DIRECT ACCESS ALLOWED.' //) +C END OF ZHOPEN. + END + SUBROUTINE ZHOPNC(IUNIT,CFILE,CFORM,CSTAT,IRLEN) + IMPLICIT NONE +C + INTEGER IUNIT,IRLEN + CHARACTER*(*) CFILE,CFORM,CSTAT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE FOR SIMPLE OPEN STATEMENTS. +C +C SEE ALSO, ZHOPEN. +C +C 2) THIS VERSION IS FOR THE SUN FOR SUN FORTRAN. +C THE FILENAME IS TAKEN FROM 'CFILE'. +C +C 3) CSTAT CAN BE 'SCRATCH', 'OLD', 'NEW', OR 'UNKNOWN'. +C CFORM CAN BE 'FORMATTED' OR 'UNFORMATTED'. +C IRLEN CAN BE ZERO (FOR SEQUENTIAL ACCESS), OR NON-ZERO (FOR DIRECT +C ACCESS INDICATING RECORD LENGTH IN TERMS OF REAL VARIABLES). +C IF IRLEN IS NEGATIVE, THE OUTPUT WILL BE IN IEEE BINARY, IF THAT +C CAPABILITY EXISTS USING STANDARD FORTRAN I/O. THIS CAPABILITY +C IS PRIMARILY TARGETED TO CRAYS, ON OTHER MACHINES -LEN AND LEN +C ARE LIKELY TO DO THE SAME THING. +C +C ON THE SUN, LEN AND -LEN BOTH GIVE IEEE FILES. +C +C 4) FOR F90 COMPILERS, DELIM='QUOTE' IS INCLUDED IN THE OPEN +C STATEMENT WHERE APPROPRIATE. +C ADDITIONALLY, FOR F90 COMPILERS: +C STATUS='NEW' IMPLIES ACTION='WRITE' +C STATUS='OLD' IMPLIES ACTION='READ' +C STATUS='SCRATCH' IMPLIES ACTION='READWRITE' +C +C 5) ALAN J. WALLCRAFT, DECEMBER 1991 AND AUGUST 1993. +C* +C********** +C + INTEGER LEN_TRIM + INTEGER IOS,NRECL +C +C OPEN FILE. +C + IF (IRLEN.EQ.0) THEN +C +C SEQUENTIAL. +C + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + IOSTAT=IOS) + ELSE +C +C UNFORMATTED DIRECT ACCESS. +C + IF (CFORM.NE.'UNFORMATTED' .AND. + + CFORM.NE.'unformatted' ) THEN + WRITE(6,9100) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF +C + IF (IRLEN.LT.0) THEN +C +C IEEE I/O. +C + NRECL = -4*IRLEN + ELSE + NRECL = 4*IRLEN + ENDIF + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ENDIF + IF (IOS.NE.0) THEN + WRITE(6,9000) IUNIT,CFILE(1:LEN_TRIM(CFILE)) + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + RETURN +C + 9000 FORMAT(// 10X,'ERROR IN ZHOPNC - CAN''T OPEN UNIT',I3,'.' / + + 10X,'CFILE = ',A //) + 9100 FORMAT(// 10X,'ERROR IN ZHOPNC (UNIT',I3.2,') -' / + + 20X,'ONLY UNFORMATTED DIRECT ACCESS ALLOWED.' //) +C END OF ZHOPNC. + END + SUBROUTINE ZHDATE(CDATE) + IMPLICIT NONE +C + CHARACTER*9 CDATE +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT RETURNS THE DATE IN 'CDATE'. +C +C 2) THE FORMAT OF CDATE NEED NOT BE IDENTICAL ON ALL MACHINES, +C BUT IT SHOULD LOOK SOMETHING LIKE, FOR EXAMPLE, '16-SEP-84', +C OR '16-Sep-84', OR ' 16/09/84' FOR THE 16TH OF SEPTEMBER 1984. +C +C 3) THIS VERSION FOR SUN WORKSTATIONS, IT USES THE +C FDATE SYSTEM ROUTINE. +C +C FDATE AS A SUN UNIX SYSTEM SUBROUTINE RETURNS A 24 +C CHARACTER STRING, FOR EAXMPLE 'Sun Sep 16 01:03:52 1984' +C FOR SOMETIME ON THE 16TH OF SEPTEMBER 1984. +C +C 4) ALAN J. WALLCRAFT, JANUARY 1989. +C* +C********** +C + CHARACTER*24 CFDATE +C + CALL FDATE(CFDATE) + CDATE = CFDATE(9:10) // '-' // CFDATE(5:7) // '-' // CFDATE(23:24) + RETURN +C END OF ZHDATE. + END + SUBROUTINE ZHSEC(SEC) + IMPLICIT NONE +C + REAL*8 SEC +C +C********** +C* +C 1) MACHINE SPECIFC ROUTINE FOR TOTAL CPU TIME UP TO THIS POINT. +C +C 2) THIS VERSION FOR THE SUN. +C +C 3) ALAN J. WALLCRAFT, PLANNING SYSTEMS INC., OCTOBER 1993. +C* +C********** +C + REAL*4 TARRAY(2) + REAL*4 ETIME +C +!MHRI SOBROUTINE NOT USED AND ETIME NOT DEFINED!!! SEC = ETIME(TARRAY) + call Cpu_Time(SEC) ! cpu_time is fortran standard + RETURN +C END OF ZHSEC. + END + INTEGER FUNCTION LEN_TRIM(CSTR) + IMPLICIT NONE +C + CHARACTER*(*) CSTR +C +C THIS FUNCTION DETERMINES THE RIGHT MOST NON-BLANK CHARACTER +C POSITION IN A STRING.. +C +C CSTR - STRING WHOSE NON-BLANK LENGTH IS TO BE DETERMINED. +C +C LEN_TRIM IS AN ELEMENTAL FUNCTION IN FORTRAN 90, +C + INTEGER I + INTEGER LEN +C + DO 110 I= LEN(CSTR), 1, -1 + IF (CSTR(I:I).NE.' ') THEN + LEN_TRIM = I + GOTO 1110 + ENDIF + 110 CONTINUE + LEN_TRIM = 0 + 1110 CONTINUE + RETURN +C END OF LEN_TRIM. + END diff --git a/sorc/hycom_tools_for_mom6.fd/zh_t3e.f b/sorc/hycom_tools_for_mom6.fd/zh_t3e.f new file mode 100644 index 00000000..a1c0f805 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/zh_t3e.f @@ -0,0 +1,287 @@ +C +C----------------------------------------------------------------------- +C +C MACHINE DEPENDENT ROUTINES. +C +C THESE ARE FOR A SINGLE T3E PROCESSOR UNDER F90 (REAL*4). +C +C----------------------------------------------------------------------- +C + SUBROUTINE ZHFLSH(IUNIT) + IMPLICIT NONE +C + INTEGER IUNIT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT FLUSHES THE OUTPUT BUFFERS OF +C LOGICAL UNIT 'IUNIT'. +C +C 2) USE ZAIOFL TO FLUSH ARRAY I/O. +C +C 3) THIS VERSION IS FOR T3E UNDER F90 (REAL*4). +C IT USES THE 'FLUSH' FORTRAN SYSTEM ROUTINE. +C +C 4) ALAN J. WALLCRAFT, SEPTEMBER 1989. +C* +C********** +C + INTEGER IOS,IRLEN + CHARACTER*240 CACC,CFILE,CFORM +C + CALL FLUSH(IUNIT,IOS) +C + IF (IOS.EQ.-1) THEN +C +C IF FLUSH DID NOT WORK, CLOSE AND RE-OPEN THE FILE. +C + INQUIRE(UNIT=IUNIT, NAME=CFILE, FORM=CFORM, + + ACCESS=CACC, RECL=IRLEN) + IF (CACC.EQ.'SEQUENTIAL') THEN + CLOSE(UNIT=IUNIT, STATUS='KEEP') + IF (CFORM.NE.'FORMATTED') THEN + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + ACTION='WRITE', POSITION='APPEND') + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + DELIM='QUOTE', RECL=4096, + + ACTION='WRITE', POSITION='APPEND') + ENDIF + ELSEIF (CACC.EQ.'DIRECT') THEN + CLOSE(UNIT=IUNIT, STATUS='KEEP') + OPEN( UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + ACTION='WRITE', ACCESS='DIRECT', RECL=IRLEN) + ENDIF + ENDIF + RETURN +C END OF ZHFLSH. + END + SUBROUTINE ZHOPEN(IUNIT,CFORM,CSTAT,IRLEN) + IMPLICIT NONE +C + INTEGER IUNIT,IRLEN + CHARACTER*(*) CFORM,CSTAT +C + INTEGER JPR + COMMON/NPROCS/ JPR + SAVE /NPROCS/ +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE FOR SIMPLE OPEN STATEMENTS. +C +C SEE ALSO, ZHOPNC. +C +C 2) THIS VERSION IS FOR THE T3E UNDER F90 (REAL*4). +C THE FILENAME IS TAKEN FROM THE ENVIRONMENT VARIABLE FOR0xx, +C WHERE xx = IUNIT, WITH DEFAULT fort.xx. +C +C 3) CSTAT CAN BE 'SCRATCH', 'OLD', 'NEW', OR 'UNKNOWN'. +C CFORM CAN BE 'FORMATTED' OR 'UNFORMATTED'. +C IRLEN CAN BE ZERO (FOR SEQUENTIAL ACCESS), OR NON-ZERO (FOR DIRECT +C ACCESS INDICATING RECORD LENGTH IN TERMS OF REAL VARIABLES). +C IF IRLEN IS NEGATIVE, THE OUTPUT WILL BE IN IEEE BINARY, IF THAT +C CAPABILITY EXISTS USING STANDARD FORTRAN I/O. THIS CAPABILITY +C IS PRIMARILY TARGETED TO CRAYS, ON OTHER MACHINES -LEN AND LEN +C ARE LIKELY TO DO THE SAME THING. +C +C ON THE T3E, LEN AND -LEN BOTH GIVE IEEE FILES. +C +C 4) FOR F90 COMPILERS, DELIM='QUOTE' IS INCLUDED IN THE OPEN +C STATEMENT WHERE APPROPRIATE. NOTE THAT THE FOLLOWING CALL +C CALL ZHOPEN(6,'FORMATTED','UNKNOWN',0) +C IS LEGAL, AND WOULD HAVE THE EFFECT OF SETTING DELIM='QUOTE' +C FOR STDOUT. IUNIT=6 IS TYPICALLY TREATED AS A SPECIAL CASE. +C ADDITIONALLY, FOR F90 COMPILERS: +C STATUS='NEW' IMPLIES ACTION='WRITE' +C STATUS='OLD' IMPLIES ACTION='READ' +C STATUS='SCRATCH' IMPLIES ACTION='READWRITE' +C +C 5) JPR FROM COMMON /NPROCS/ IS THE EXPECTED JPR PARAMETER IN THE +C TARGET OCEAN MODEL. +C +C 6) ALAN J. WALLCRAFT, NRL, MAY 1997. +C* +C********** +C + INTEGER IOS,NRECL + CHARACTER CFILE*240,CENV*6,CACT*9 +C + INTEGER IBLOCK,ICACHE,JCACHE + CHARACTER CASN*40 +C +C GET FILENAME. +C + WRITE(CENV,1000) IUNIT + CFILE = ' ' + CALL GETENV(CENV,CFILE) + IF (CFILE.EQ.' ') THEN + WRITE(CFILE,1100) IUNIT + ENDIF +C +C OPEN FILE. +C + IF (CSTAT.EQ.'OLD' .OR. + + CSTAT.EQ.'old' ) THEN + CACT = 'READ' + ELSEIF (CSTAT.EQ.'NEW' .OR. + + CSTAT.EQ.'new' ) THEN + CACT = 'WRITE' + ELSE + CACT = 'READWRITE' + ENDIF +C + IF (IRLEN.EQ.0) THEN +C +C SEQUENTIAL (UNFORMATTED USES IEEE I/O WITH F77/UNIX BLOCKING). +C + IF (CFORM.EQ.'UNFORMATTED' .OR. + + CFORM.EQ.'unformatted' ) THEN + CALL ASNUNIT(IUNIT,'-F f77',IOS) + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACTION=CACT, IOSTAT=IOS) + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, IOSTAT=IOS) + ENDIF + ELSE + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACTION=CACT, DELIM='QUOTE', RECL=4096, IOSTAT=IOS) + ELSEIF (IUNIT.NE.6) THEN + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, DELIM='QUOTE', RECL=4096, IOSTAT=IOS) + ELSE + IOS = 0 +* OPEN(UNIT=6, DELIM='QUOTE', RECL=4096, IOSTAT=IOS) + ENDIF + ENDIF + ELSE +C +C UNFORMATTED DIRECT ACCESS. +C + IF (CFORM.NE.'UNFORMATTED' .AND. + + CFORM.NE.'unformatted' ) THEN + WRITE(6,9100) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF +C + IF (IRLEN.LT.0) THEN +C +C IEEE I/O. +C + NRECL = -4*IRLEN + ELSE + NRECL = 4*IRLEN + ENDIF + IF (MOD(NRECL,4096).EQ.0) THEN + ICACHE = NRECL/4096 + IF (MOD(ICACHE,JPR).EQ.0) THEN + JCACHE = ICACHE/JPR + IF (MOD(JCACHE,3).EQ.0) THEN + IBLOCK = JCACHE/3 + ELSEIF (MOD(JCACHE,2).EQ.0) THEN + IBLOCK = JCACHE/2 + ELSE + IBLOCK = JCACHE + ENDIF + WRITE(CASN,8000) ICACHE,IBLOCK + 8000 FORMAT('-F cachea:',I4.4,':1:0 -p6-63 -q',I4.4) + CALL ASNUNIT(IUNIT,CASN,IOS) + ENDIF + ENDIF + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACTION=CACT, ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ENDIF + ENDIF + IF (IOS.NE.0) THEN + WRITE(6,9000) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF + RETURN +C + 1000 FORMAT('FOR',I3.3) + 1100 FORMAT('fort.',I2.2) + 9000 FORMAT(// 10X,'ERROR IN ZHOPEN - CAN''T OPEN UNIT',I3,'.' //) + 9100 FORMAT(// 10X,'ERROR IN ZHOPEN (UNIT',I3.2,') -' / + + 20X,'ONLY UNFORMATTED DIRECT ACCESS ALLOWED.' //) +C END OF ZHOPEN. + END + SUBROUTINE ZHDATE(CDATE) + IMPLICIT NONE +C + CHARACTER*9 CDATE +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT RETURNS THE DATE IN 'CDATE'. +C +C 2) THE FORMAT OF CDATE NEED NOT BE IDENTICAL ON ALL MACHINES, +C BUT IT SHOULD LOOK SOMETHING LIKE, FOR EXAMPLE, '16-SEP-84', +C OR '16-Sep-84', OR ' 16/09/84' FOR THE 16TH OF SEPTEMBER 1984. +C +C 3) THIS VERSION FOR THE T3E, IT USES THE +C DATE SYSTEM ROUTINE. +C +C DATE AS A CRAY SYSTEM SUBROUTINE RETURNS, FOR EXAMPLE +C '16/09/84' FOR THE 16TH OF SEPTEMBER 1984. +C +C 4) ALAN J. WALLCRAFT, JANUARY 1989. +C* +C********** +C + CDATE(1:1) = ' ' + CALL DATE(CDATE(2:9)) + RETURN +C END OF ZHDATE. + END + SUBROUTINE ZHSEC(SEC) + IMPLICIT NONE +C + REAL*8 SEC +C +C********** +C* +C 1) MACHINE SPECIFC ROUTINE FOR TOTAL CPU TIME UP TO THIS POINT. +C +C 2) THIS VERSION FOR THE T3E UNDER F90 (REAL*4). +C +C 3) ALAN J. WALLCRAFT, OCTOBER 1993. +C* +C********** +C + INTEGER IRTC +C + SEC = IRTC() * 3.333333333E-9 + RETURN +C END OF ZHSEC. + END + SUBROUTINE GETENV(CNAME, CVALUE) + IMPLICIT NONE +C + CHARACTER*(*) CNAME,CVALUE +C +C THIS SUBROUTINE PROVIDES GETENV FUNCTIONALITY +C ON THE T3E, USING PXFGETENV. +C + INTEGER INAME,IVALUE,IERR +C + INAME = 0 + IERR = 0 + CALL PXFGETENV(CNAME,INAME, CVALUE,IVALUE, IERR) + IF (IERR.NE.0) THEN + CVALUE = ' ' + ENDIF + RETURN +C END OF GETENV. + END diff --git a/sorc/hycom_tools_for_mom6.fd/zh_x1.f b/sorc/hycom_tools_for_mom6.fd/zh_x1.f new file mode 100644 index 00000000..1a8ac4d4 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/zh_x1.f @@ -0,0 +1,322 @@ +C +C----------------------------------------------------------------------- +C +C MACHINE DEPENDENT ROUTINES. +C +C THESE ARE FOR A CRAY_X1 UNDER FTN (REAL*4). +C +C----------------------------------------------------------------------- +C + SUBROUTINE ZHFLSH(IUNIT) + IMPLICIT NONE +C + INTEGER IUNIT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT FLUSHES THE OUTPUT BUFFERS OF +C LOGICAL UNIT 'IUNIT'. +C +C 2) USE ZAIOFL TO FLUSH ARRAY I/O. +C +C 3) THIS VERSION IS FOR CRAY_X1 UNDER FTN (REAL*4). +C IT USES THE 'FLUSH' FORTRAN SYSTEM ROUTINE. +C +C 4) ALAN J. WALLCRAFT, SEPTEMBER 1989. +C* +C********** +C + INTEGER IOS,IRLEN + CHARACTER*240 CACC,CFILE,CFORM +C + CALL FLUSH(IUNIT,IOS) +C + IF (IOS.EQ.-1) THEN +C +C IF FLUSH DID NOT WORK, CLOSE AND RE-OPEN THE FILE. +C + INQUIRE(UNIT=IUNIT, NAME=CFILE, FORM=CFORM, + + ACCESS=CACC, RECL=IRLEN) + IF (CACC.EQ.'SEQUENTIAL') THEN + CLOSE(UNIT=IUNIT, STATUS='KEEP') + IF (CFORM.NE.'FORMATTED') THEN + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + ACTION='WRITE', POSITION='APPEND') + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + DELIM='QUOTE', RECL=4096, + + ACTION='WRITE', POSITION='APPEND') + ENDIF + ELSEIF (CACC.EQ.'DIRECT') THEN + CLOSE(UNIT=IUNIT, STATUS='KEEP') + OPEN( UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + ACTION='WRITE', ACCESS='DIRECT', RECL=IRLEN) + ENDIF + ENDIF + RETURN +C END OF ZHFLSH. + END + SUBROUTINE ZHOPEN(IUNIT,CFORM,CSTAT,IRLEN) + IMPLICIT NONE +C + INTEGER IUNIT,IRLEN + CHARACTER*(*) CFORM,CSTAT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE FOR SIMPLE OPEN STATEMENTS. +C +C SEE ALSO, ZHOPNC. +C +C 2) THIS VERSION IS FOR THE CRAY_X1 UNDER FTN (REAL*4). +C THE FILENAME IS TAKEN FROM THE ENVIRONMENT VARIABLE FOR0xx, +C WHERE xx = IUNIT, WITH DEFAULT fort.xx. +C +C 3) CSTAT CAN BE 'SCRATCH', 'OLD', 'NEW', OR 'UNKNOWN'. +C CFORM CAN BE 'FORMATTED' OR 'UNFORMATTED'. +C IRLEN CAN BE ZERO (FOR SEQUENTIAL ACCESS), OR NON-ZERO (FOR DIRECT +C ACCESS INDICATING RECORD LENGTH IN TERMS OF REAL VARIABLES). +C IF IRLEN IS NEGATIVE, THE OUTPUT WILL BE IN IEEE BINARY, IF THAT +C CAPABILITY EXISTS USING STANDARD FORTRAN I/O. THIS CAPABILITY +C IS PRIMARILY TARGETED TO CRAYS, ON OTHER MACHINES -LEN AND LEN +C ARE LIKELY TO DO THE SAME THING. +C +C ON THE CRAY_X1, LEN AND -LEN BOTH GIVE IEEE FILES. +C +C 4) FOR F90 COMPILERS, DELIM='QUOTE' IS INCLUDED IN THE OPEN +C STATEMENT WHERE APPROPRIATE. NOTE THAT THE FOLLOWING CALL +C CALL ZHOPEN(6,'FORMATTED','UNKNOWN',0) +C IS LEGAL, AND WOULD HAVE THE EFFECT OF SETTING DELIM='QUOTE' +C FOR STDOUT. IUNIT=6 IS TYPICALLY TREATED AS A SPECIAL CASE. +C ADDITIONALLY, FOR F90 COMPILERS: +C STATUS='NEW' IMPLIES ACTION='WRITE' +C STATUS='OLD' IMPLIES ACTION='READ' +C STATUS='SCRATCH' IMPLIES ACTION='READWRITE' +C +C 5) ALAN J. WALLCRAFT, DECEMBER 1991 AND AUGUST 1993. +C* +C********** +C + INTEGER IOS,NRECL + CHARACTER CFILE*240,CENV*6 +C +C GET FILENAME. +C + WRITE(CENV,1000) IUNIT + CFILE = ' ' + CALL GETENV(CENV,CFILE) + IF (CFILE.EQ.' ') THEN + WRITE(CFILE,1100) IUNIT + ENDIF +C +C OPEN FILE. +C + IF (IRLEN.EQ.0) THEN +C +C SEQUENTIAL. +C + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + IOSTAT=IOS) + ELSEIF (IUNIT.NE.6) THEN + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + IOSTAT=IOS) + ELSE + IOS = 0 + ENDIF + ELSE +C +C UNFORMATTED DIRECT ACCESS. +C + IF (CFORM.NE.'UNFORMATTED' .AND. + + CFORM.NE.'unformatted' ) THEN + WRITE(6,9100) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF +C + IF (IRLEN.LT.0) THEN +C +C IEEE I/O. +C + NRECL = -4*IRLEN + ELSE + NRECL = 4*IRLEN + ENDIF + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ENDIF + ENDIF + IF (IOS.NE.0) THEN + WRITE(6,9000) IUNIT + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + RETURN +C + 1000 FORMAT('FOR',I3.3) + 1100 FORMAT('fort.',I2.2) + 9000 FORMAT(/ / 10X,'ERROR IN ZHOPEN - CAN''T OPEN UNIT',I3,'.' / /) + 9100 FORMAT(/ / 10X,'ERROR IN ZHOPEN (UNIT',I3.2,') -' / + + 20X,'ONLY UNFORMATTED DIRECT ACCESS ALLOWED.' / /) +C END OF ZHOPEN. + END + SUBROUTINE ZHOPNC(IUNIT,CFILE,CFORM,CSTAT,IRLEN) + IMPLICIT NONE +C + INTEGER IUNIT,IRLEN + CHARACTER*(*) CFILE,CFORM,CSTAT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE FOR SIMPLE OPEN STATEMENTS. +C +C SEE ALSO, ZHOPEN. +C +C 2) THIS VERSION IS FOR THE CRAY_X1 UNDER FTN (REAL*4). +C THE FILENAME IS TAKEN FROM 'CFILE'. +C +C 3) CSTAT CAN BE 'SCRATCH', 'OLD', 'NEW', OR 'UNKNOWN'. +C CFORM CAN BE 'FORMATTED' OR 'UNFORMATTED'. +C IRLEN CAN BE ZERO (FOR SEQUENTIAL ACCESS), OR NON-ZERO (FOR DIRECT +C ACCESS INDICATING RECORD LENGTH IN TERMS OF REAL VARIABLES). +C IF IRLEN IS NEGATIVE, THE OUTPUT WILL BE IN IEEE BINARY, IF THAT +C CAPABILITY EXISTS USING STANDARD FORTRAN I/O. THIS CAPABILITY +C IS PRIMARILY TARGETED TO CRAYS, ON OTHER MACHINES -LEN AND LEN +C ARE LIKELY TO DO THE SAME THING. +C +C ON THE CRAY_X1, LEN AND -LEN BOTH GIVE IEEE FILES. +C +C 4) FOR F90 COMPILERS, DELIM='QUOTE' IS INCLUDED IN THE OPEN +C STATEMENT WHERE APPROPRIATE. +C ADDITIONALLY, FOR F90 COMPILERS: +C STATUS='NEW' IMPLIES ACTION='WRITE' +C STATUS='OLD' IMPLIES ACTION='READ' +C STATUS='SCRATCH' IMPLIES ACTION='READWRITE' +C +C 5) ALAN J. WALLCRAFT, DECEMBER 1991 AND AUGUST 1993. +C* +C********** +C + INTEGER LEN_TRIM + INTEGER IOS,NRECL +C +C OPEN FILE. +C + IF (IRLEN.EQ.0) THEN +C +C SEQUENTIAL. +C + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + IOSTAT=IOS) + ELSE +C +C UNFORMATTED DIRECT ACCESS. +C + IF (CFORM.NE.'UNFORMATTED' .AND. + + CFORM.NE.'unformatted' ) THEN + WRITE(6,9100) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF +C + IF (IRLEN.LT.0) THEN +C +C IEEE I/O. +C + NRECL = -4*IRLEN + ELSE + NRECL = 4*IRLEN + ENDIF + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ENDIF + IF (IOS.NE.0) THEN + WRITE(6,9000) IUNIT,CFILE(1:LEN_TRIM(CFILE)) + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + RETURN +C + 9000 FORMAT(/ / 10X,'ERROR IN ZHOPNC - CAN''T OPEN UNIT',I3,'.' / + + 10X,'CFILE = ',A / /) + 9100 FORMAT(/ / 10X,'ERROR IN ZHOPNC (UNIT',I3.2,') -' / + + 20X,'ONLY UNFORMATTED DIRECT ACCESS ALLOWED.' / /) +C END OF ZHOPNC. + END + SUBROUTINE ZHDATE(CDATE) + IMPLICIT NONE +C + CHARACTER*9 CDATE +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT RETURNS THE DATE IN 'CDATE'. +C +C 2) THE FORMAT OF CDATE NEED NOT BE IDENTICAL ON ALL MACHINES, +C BUT IT SHOULD LOOK SOMETHING LIKE, FOR EXAMPLE, '16-SEP-84', +C OR '16-Sep-84', OR ' 16/09/84' FOR THE 16TH OF SEPTEMBER 1984. +C +C 3) THIS VERSION FOR THE CRAY_X1, IT USES THE +C DATE SYSTEM ROUTINE. +C +C DATE AS A CRAY SYSTEM SUBROUTINE RETURNS, FOR EXAMPLE +C '16/09/84' FOR THE 16TH OF SEPTEMBER 1984. +C +C 4) ALAN J. WALLCRAFT, JANUARY 1989. +C* +C********** +C + CDATE(1:1) = ' ' + CALL DATE(CDATE(2:9)) + RETURN +C END OF ZHDATE. + END + SUBROUTINE ZHSEC(SEC) + IMPLICIT NONE +C + REAL*8 SEC +C +C********** +C* +C 1) MACHINE SPECIFC ROUTINE FOR TOTAL CPU TIME UP TO THIS POINT. +C +C 2) THIS VERSION FOR THE CRAY_X1 UNDER FTN (REAL*4). +C +C 3) ALAN J. WALLCRAFT, OCTOBER 1993. +C* +C********** +C + REAL*4 SECOND +C + SEC = SECOND() + RETURN +C END OF ZHSEC. + END + SUBROUTINE GETENV(CNAME, CVALUE) + IMPLICIT NONE +C + CHARACTER*(*) CNAME,CVALUE +C +C THIS SUBROUTINE PROVIDES GETENV FUNCTIONALITY +C ON THE CRAY_X1 USING PXFGETENV. +C + INTEGER INAME,IVALUE,IERR +C + INAME = 0 + IERR = 0 + CALL PXFGETENV(CNAME,INAME, CVALUE,IVALUE, IERR) + IF (IERR.NE.0) THEN + CVALUE = ' ' + ENDIF + RETURN +C END OF GETENV. + END diff --git a/sorc/hycom_tools_for_mom6.fd/zh_ymp.f b/sorc/hycom_tools_for_mom6.fd/zh_ymp.f new file mode 100644 index 00000000..25447496 --- /dev/null +++ b/sorc/hycom_tools_for_mom6.fd/zh_ymp.f @@ -0,0 +1,419 @@ +C +C----------------------------------------------------------------------- +C +C MACHINE DEPENDENT ROUTINES. +C +C THESE ARE FOR THE CRAY UNDER UNICOS 8.0. +C +C----------------------------------------------------------------------- +C + SUBROUTINE ZHFLSH(IUNIT) + IMPLICIT NONE +C + INTEGER IUNIT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT FLUSHES THE OUTPUT BUFFERS OF +C LOGICAL UNIT 'IUNIT'. +C +C 2) USE ZAIOFL TO FLUSH ARRAY I/O. +C +C 3) THIS VERSION IS FOR CRAY C90 UNDER F90. +C IT USES THE 'FLUSH' FORTRAN SYSTEM ROUTINE. +C +C 4) ALAN J. WALLCRAFT, SEPTEMBER 1989. +C* +C********** +C + INTEGER IOS,IRLEN + CHARACTER*240 CACC,CFILE,CFORM +C + CALL FLUSH(IUNIT,IOS) +C + IF (IOS.EQ.-1) THEN +C +C IF FLUSH DID NOT WORK, CLOSE AND RE-OPEN THE FILE. +C + INQUIRE(UNIT=IUNIT, NAME=CFILE, FORM=CFORM, + + ACCESS=CACC, RECL=IRLEN) + CLOSE( UNIT=IUNIT, STATUS='KEEP') + IF (CACC.NE.'DIRECT') THEN + IF (CFORM.NE.'FORMATTED') THEN + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + ACTION='WRITE', POSITION='APPEND') + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + DELIM='QUOTE', RECL=4096, + + ACTION='WRITE', POSITION='APPEND') + ENDIF + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS='OLD', + + ACTION='WRITE', ACCESS='DIRECT', RECL=IRLEN) + ENDIF + ENDIF + RETURN +C END OF ZHFLSH. + END + SUBROUTINE ZHOPEN(IUNIT,CFORM,CSTAT,IRLEN) + IMPLICIT NONE +C + INTEGER IUNIT,IRLEN + CHARACTER*(*) CFORM,CSTAT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE FOR SIMPLE OPEN STATEMENTS. +C +C SEE ALSO, ZHOPNC. +C +C 2) THIS VERSION IS FOR THE CRAY C90 UNDER F90. +C THE FILENAME IS TAKEN FROM THE ENVIRONMENT VARIABLE FOR0xx, +C WHERE xx = IUNIT, WITH DEFAULT fort.xx. +C +C 3) CSTAT CAN BE 'SCRATCH', 'OLD', 'NEW', OR 'UNKNOWN'. +C CFORM CAN BE 'FORMATTED' OR 'UNFORMATTED'. +C IRLEN CAN BE ZERO (FOR SEQUENTIAL ACCESS), OR NON-ZERO (FOR DIRECT +C ACCESS INDICATING RECORD LENGTH IN TERMS OF REAL VARIABLES). +C IF IRLEN IS NEGATIVE, THE OUTPUT WILL BE IN IEEE BINARY, IF THAT +C CAPABILITY EXISTS USING STANDARD FORTRAN I/O. THIS CAPABILITY +C IS PRIMARILY TARGETED TO CRAYS, ON OTHER MACHINES -LEN AND LEN +C ARE LIKELY TO DO THE SAME THING. +C +C ON THE SUN, LEN AND -LEN BOTH GIVE IEEE FILES. +C +C 4) FOR F90 COMPILERS, DELIM='QUOTE' IS INCLUDED IN THE OPEN +C STATEMENT WHERE APPROPRIATE. NOTE THAT THE FOLLOWING CALL +C CALL ZHOPEN(6,'FORMATTED','UNKNOWN',0) +C IS LEGAL, AND WOULD HAVE THE EFFECT OF SETTING DELIM='QUOTE' +C FOR STDOUT. IUNIT=6 IS TYPICALLY TREATED AS A SPECIAL CASE. +C ADDITIONALLY, FOR F90 COMPILERS: +C STATUS='NEW' IMPLIES ACTION='WRITE' +C STATUS='OLD' IMPLIES ACTION='READ' +C STATUS='SCRATCH' IMPLIES ACTION='READWRITE' +C +C 5) ALAN J. WALLCRAFT, DECEMBER 1991 AND AUGUST 1993. +C* +C********** +C + INTEGER IOS,NRECL + CHARACTER CFILE*240,CENV*6,CACT*9,CFCONV*8,CFSPEC*18 +C +C GET FILENAME. +C + WRITE(CENV,1000) IUNIT + CFILE = ' ' + CALL GETENV(CENV,CFILE) + IF (CFILE.EQ.' ') THEN + WRITE(CFILE,1100) IUNIT + ENDIF +C +C OPEN FILE. +C + IF (CSTAT.EQ.'OLD' .OR. + + CSTAT.EQ.'old' ) THEN + CACT = 'READ' + ELSEIF (CSTAT.EQ.'NEW' .OR. + + CSTAT.EQ.'new' ) THEN + CACT = 'WRITE' + ELSE + CACT = 'READWRITE' + ENDIF +C + IF (IRLEN.EQ.0) THEN +C +C SEQUENTIAL (UNFORMATTED USES IEEE I/O WITH F77/UNIX BLOCKING). +C + IF (CFORM.EQ.'UNFORMATTED' .OR. + + CFORM.EQ.'unformatted' ) THEN + CALL ASNUNIT(IUNIT,'-F f77 -N ieee',IOS) + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACTION=CACT, IOSTAT=IOS) + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, IOSTAT=IOS) + ENDIF + ELSE + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACTION=CACT, DELIM='QUOTE', RECL=4096, IOSTAT=IOS) + ELSEIF (IUNIT.NE.6) THEN + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, DELIM='QUOTE', RECL=4096, IOSTAT=IOS) + ELSE + OPEN(UNIT=6, DELIM='QUOTE', RECL=4096, IOSTAT=IOS) + ENDIF + ENDIF + ELSE +C +C UNFORMATTED DIRECT ACCESS. +C + IF (CFORM.NE.'UNFORMATTED' .AND. + + CFORM.NE.'unformatted' ) THEN + WRITE(6,9100) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF +C + IF (IRLEN.LT.0) THEN +C +C IEEE I/O. +C + CFCONV = '-N ieee' + NRECL = -4*IRLEN + ELSE + CFCONV = ' ' + NRECL = 8*IRLEN + ENDIF + IF (MOD(NRECL,16384).EQ.0 .AND. NRECL.GT.16384*4) THEN +C +C NO BUFFERS FOR LARGE, WELL-FORMED RECORDS. +C + CFSPEC = '-F syscall' + ELSE +C +C CACHEA FOR SMALL RECORDS, +C 16 PAGES OF 8 BLOCKS WITH READ-AHEAD OF TWO. +C + CFSPEC = '-F cachea:8:16:2' + ENDIF + CALL ASNUNIT(IUNIT,CFSPEC//CFCONV,IOS) + IF (IOS.NE.0) THEN + WRITE(6,9050) IUNIT + WRITE(6,*) 'CFSPEC = ',CFSPEC + WRITE(6,*) 'CFCONV = ',CFCONV + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + IF (CSTAT.EQ.'SCRATCH' .OR. + + CSTAT.EQ.'scratch' ) THEN + OPEN(UNIT=IUNIT, FORM=CFORM, STATUS='SCRATCH', + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ENDIF + ENDIF + IF (IOS.NE.0) THEN + WRITE(6,9000) IUNIT + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + RETURN +C + 1000 FORMAT('FOR',I3.3) + 1100 FORMAT('fort.',I2.2) + 9000 FORMAT(// 10X,'ERROR IN ZHOPEN - CAN''T OPEN UNIT',I3,'.' //) + 9050 FORMAT(// 10X,'ERROR IN ZHOPEN - CAN''T ASNUNIT',I3,'.' //) + 9100 FORMAT(// 10X,'ERROR IN ZHOPEN (UNIT',I3.2,') -' / + + 20X,'ONLY UNFORMATTED DIRECT ACCESS ALLOWED.' //) +C END OF ZHOPEN. + END + SUBROUTINE ZHOPNC(IUNIT,CFILE,CFORM,CSTAT,IRLEN) + IMPLICIT NONE +C + INTEGER IUNIT,IRLEN + CHARACTER*(*) CFILE,CFORM,CSTAT +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE FOR SIMPLE OPEN STATEMENTS. +C +C SEE ALSO, ZHOPEN. +C +C 2) THIS VERSION IS FOR THE CRAY C90 UNDER F90. +C THE FILENAME IS TAKEN FROM 'CFILE'. +C +C 3) CSTAT CAN BE 'SCRATCH', 'OLD', 'NEW', OR 'UNKNOWN'. +C CFORM CAN BE 'FORMATTED' OR 'UNFORMATTED'. +C IRLEN CAN BE ZERO (FOR SEQUENTIAL ACCESS), OR NON-ZERO (FOR DIRECT +C ACCESS INDICATING RECORD LENGTH IN TERMS OF REAL VARIABLES). +C IF IRLEN IS NEGATIVE, THE OUTPUT WILL BE IN IEEE BINARY, IF THAT +C CAPABILITY EXISTS USING STANDARD FORTRAN I/O. THIS CAPABILITY +C IS PRIMARILY TARGETED TO CRAYS, ON OTHER MACHINES -LEN AND LEN +C ARE LIKELY TO DO THE SAME THING. +C +C ON THE SUN, LEN AND -LEN BOTH GIVE IEEE FILES. +C +C 4) FOR F90 COMPILERS, DELIM='QUOTE' IS INCLUDED IN THE OPEN +C STATEMENT WHERE APPROPRIATE. +C ADDITIONALLY, FOR F90 COMPILERS: +C STATUS='NEW' IMPLIES ACTION='WRITE' +C STATUS='OLD' IMPLIES ACTION='READ' +C STATUS='SCRATCH' IMPLIES ACTION='READWRITE' +C +C 5) ALAN J. WALLCRAFT, DECEMBER 1991 AND AUGUST 1993. +C* +C********** +C + INTEGER LEN_TRIM + INTEGER IOS,NRECL + CHARACTER CACT*9,CFCONV*8,CFSPEC*18 +C +C OPEN FILE. +C + IF (CSTAT.EQ.'OLD' .OR. + + CSTAT.EQ.'old' ) THEN + CACT = 'READ' + ELSEIF (CSTAT.EQ.'NEW' .OR. + + CSTAT.EQ.'new' ) THEN + CACT = 'WRITE' + ELSE + CACT = 'READWRITE' + ENDIF +C + IF (IRLEN.EQ.0) THEN +C +C SEQUENTIAL (UNFORMATTED USES IEEE I/O WITH F77/UNIX BLOCKING). +C + IF (CFORM.EQ.'UNFORMATTED' .OR. + + CFORM.EQ.'unformatted' ) THEN + CALL ASNUNIT(IUNIT,'-F f77 -N ieee',IOS) + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, IOSTAT=IOS) + ELSE + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, DELIM='QUOTE', RECL=4096, IOSTAT=IOS) + ENDIF + ELSE +C +C UNFORMATTED DIRECT ACCESS. +C + IF (CFORM.NE.'UNFORMATTED' .AND. + + CFORM.NE.'unformatted' ) THEN + WRITE(6,9100) IUNIT + CALL ZHFLSH(6) + STOP + ENDIF +C + IF (IRLEN.LT.0) THEN +C +C IEEE I/O. +C + CFCONV = '-N ieee' + NRECL = -4*IRLEN + ELSE + CFCONV = ' ' + NRECL = 8*IRLEN + ENDIF + IF (MOD(NRECL,16384).EQ.0.AND.NRECL.GT.16384*4) THEN +C +C NO BUFFERS FOR LARGE, WELL-FORMED RECORDS. +C + CFSPEC = '-F syscall' + ELSE +C +C CACHEA FOR SMALL RECORDS, +C 16 PAGES OF 8 BLOCKS WITH READ-AHEAD OF TWO. +C + CFSPEC = '-F cachea:8:16:2' + ENDIF + CALL ASNUNIT(IUNIT,CFSPEC//CFCONV,IOS) + IF (IOS.NE.0) THEN + WRITE(6,9050) IUNIT + WRITE(6,*) 'CFSPEC = ',CFSPEC + WRITE(6,*) 'CFCONV = ',CFCONV + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + OPEN(UNIT=IUNIT, FILE=CFILE, FORM=CFORM, STATUS=CSTAT, + + ACTION=CACT, ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) + ENDIF + IF (IOS.NE.0) THEN + WRITE(6,9000) IUNIT,CFILE(1:LEN_TRIM(CFILE)) + WRITE(6,*) 'IOS = ',IOS + CALL ZHFLSH(6) + STOP + ENDIF + RETURN +C + 1000 FORMAT('FOR',I3.3) + 1100 FORMAT('fort.',I2.2) + 9000 FORMAT(// 10X,'ERROR IN ZHOPNC - CAN''T OPEN UNIT',I3,'.' / + + 10X,'CFILE = ',A //) + 9050 FORMAT(// 10X,'ERROR IN ZHOPEN - CAN''T ASNUNIT',I3,'.' //) + 9100 FORMAT(// 10X,'ERROR IN ZHOPNC (UNIT',I3.2,') -' / + + 20X,'ONLY UNFORMATTED DIRECT ACCESS ALLOWED.' //) +C END OF ZHOPNC. + END + SUBROUTINE ZHDATE(CDATE) + IMPLICIT NONE +C + CHARACTER*9 CDATE +C +C********** +C* +C 1) MACHINE SPECIFIC ROUTINE THAT RETURNS THE DATE IN 'CDATE'. +C +C 2) THE FORMAT OF CDATE NEED NOT BE IDENTICAL ON ALL MACHINES, +C BUT IT SHOULD LOOK SOMETHING LIKE, FOR EXAMPLE, '16-SEP-84', +C OR '16-Sep-84', OR ' 16/09/84' FOR THE 16TH OF SEPTEMBER 1984. +C +C 3) THIS VERSION FOR THE CRAY C90 UNDER F90, IT USES THE +C DATE SYSTEM ROUTINE. +C +C DATE AS A CRAY SYSTEM SUBROUTINE RETURNS, FOR EXAMPLE +C '16/09/84' FOR THE 16TH OF SEPTEMBER 1984. +C +C 4) ALAN J. WALLCRAFT, JANUARY 1989. +C* +C********** +C + CDATE(1:1) = ' ' + CALL DATE(CDATE(2:9)) + RETURN +C END OF ZHDATE. + END + SUBROUTINE ZHSEC(SEC) + IMPLICIT NONE +C + REAL*8 SEC +C +C********** +C* +C 1) MACHINE SPECIFC ROUTINE FOR TOTAL CPU TIME UP TO THIS POINT. +C +C 2) THIS VERSION FOR THE SUN. +C +C 3) ALAN J. WALLCRAFT, PLANNING SYSTEMS INC., OCTOBER 1993. +C* +C********** +C + REAL*8 SECOND +C + SEC = SECOND() + RETURN +C END OF ZHSEC. + END + INTEGER FUNCTION LEN_TRIM(CSTR) + IMPLICIT NONE +C + CHARACTER*(*) CSTR +C +C THIS FUNCTION DETERMINES THE RIGHT MOST NON-BLANK CHARACTER +C POSITION IN A STRING.. +C +C CSTR - STRING WHOSE NON-BLANK LENGTH IS TO BE DETERMINED. +C +C LEN_TRIM IS AN ELEMENTAL FUNCTION IN FORTRAN 90, +C + INTEGER I + INTEGER LEN +C + DO 110 I= LEN(CSTR), 1, -1 + IF (CSTR(I:I).NE.' ') THEN + LEN_TRIM = I + GOTO 1110 + ENDIF + 110 CONTINUE + LEN_TRIM = 0 + 1110 CONTINUE + RETURN +C END OF LEN_TRIM. + END