Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_and_exec_oxf =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_and_exec_oxf (revision 1439) +++ (revision ) @@ -1,11 +1,0 @@ -#! /bin/bash - -./compile_oxf - -echo I assume input_diagfi is correctly linked -echo I assume WPSFEED folder exists or is linked - -echo 1 | create_readmeteo.exe -readmeteo.exe < readmeteo.def - -ncrcat -O -v totoptdep input_diagfi.nc input_totoptdep.nc Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_oxf =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_oxf (revision 1439) +++ (revision ) @@ -1,22 +1,0 @@ -#! /bin/bash - -\rm fix_no_info.inc 2> /dev/null -\rm *.o 2> /dev/null -\rm readmeteo_tmp.F90 2> /dev/null -cp fix_no_info_OXF.inc fix_no_info.inc - -sed s/"PARAMETER :: ident='LMD'"/"PARAMETER :: ident='OXF'"/g readmeteo.F90 > readmeteo_tmp.F90 - -pgf90 -byteswapio readmeteo_tmp.F90 \ --L$NETCDF/lib -lnetcdf \ --I$NETCDF/include \ --o readmeteo.exe - -pgf90 create_readmeteo.F90 \ --L$NETCDF/lib -lnetcdf \ --I$NETCDF/include \ --o create_readmeteo.exe - -\rm fix_no_info.inc 2> /dev/null -\rm readmeteo_tmp.F90 2> /dev/null -\rm *.o 2> /dev/null Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_oxf_g95 =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_oxf_g95 (revision 1439) +++ (revision ) @@ -1,22 +1,0 @@ -#! /bin/bash - -\rm fix_no_info.inc 2> /dev/null -\rm *.o 2> /dev/null -\rm readmeteo_tmp.F90 2> /dev/null -cp fix_no_info_OXF.inc fix_no_info.inc - -sed s/"PARAMETER :: ident='LMD'"/"PARAMETER :: ident='OXF'"/g readmeteo.F90 > readmeteo_tmp.F90 - -g95 -fno-second-underscore readmeteo_tmp.F90 \ --L$NETCDF/lib -lnetcdf \ --I$NETCDF/include \ --o readmeteo.exe - -g95 -fno-second-underscore create_readmeteo.F90 \ --L$NETCDF/lib -lnetcdf \ --I$NETCDF/include \ --o create_readmeteo.exe - -\rm fix_no_info.inc 2> /dev/null -\rm readmeteo_tmp.F90 2> /dev/null -\rm *.o 2> /dev/null Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/concatoxf.pro =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/concatoxf.pro (revision 1439) +++ (revision ) @@ -1,136 +1,0 @@ -pro concatoxf - -;------------------------------------ -file1='input_totoptdep1.nc' -file2='input_totoptdep2.nc' -;------------------------------------ - -cdfid1 = ncdf_open(file1) -cdfid2 = ncdf_open(file2) - -; coordinates -varid=ncdf_varid(cdfid1,'lat') - ncdf_varget, cdfid1, varid, lati -lat=n_elements(lati) - -varid=ncdf_varid(cdfid1,'lon') - ncdf_varget, cdfid1, varid, longi -lon=n_elements(longi) - -varid=ncdf_varid(cdfid1,'time') - ncdf_varget, cdfid1, varid, time1 -varid=ncdf_varid(cdfid2,'time') - ncdf_varget, cdfid2, varid, time2 -t1=n_elements(time1) -t2=n_elements(time2) - -; concat time coordinates -ntime=t1+t2 -timear=fltarr(ntime) -timear(0:t1-1)=time1(0:t1-1) -timear(t1:ntime-1)=time2(0:t2-1)+time1(t1-1) - - -; field -varid=ncdf_varid(cdfid1,'totoptdep') - ncdf_varget, cdfid1, varid, totoptdep1 -varid=ncdf_varid(cdfid2,'totoptdep') - ncdf_varget, cdfid2, varid, totoptdep2 - - -;loadct, 3 -;!p.multi=[0,1,3] -;ref=totoptdep1(*,*,tata) -;smo=smooth(ref,toto) -;err=100.*abs(ref-smo)/ref -;err=(ref-smo) -;print, max(err), ref(where(err eq max(err))) -;contour, err, nlevels=40, /cell_fill -;contour, ref, nlevels=40, /cell_fill -;contour, smo, nlevels=40, /cell_fill -;stop - - - -; concat field -totoptdep=fltarr(lon,lat,ntime) -totoptdep(*,*,0:t1-1)=totoptdep1(*,*,*) -totoptdep(*,*,t1:ntime-1)=totoptdep2(*,*,*) -fieldname1='totoptdep' -outfieldstring1='Total dust optical depth (tauref/700Pa)' -unit1='Opacity/Pa' - - - -;;;********************************************************************* - -; -; SMOOTH wrt SPACE -; -totoptdep=smooth(totoptdep,2) - - -; -; INTERPOLATE wrt TIME -; - -which_utc_is_TES=3 -;; for Hellas 8AM is 02PM LT, i.e. TES measurements -;; IDL convention, starts 0 !!!! - - -nday=floor(float(ntime-1)/12.)-1 - -ref=fltarr(nday) -subsarr=fltarr(nday) - -for i=0, n_elements(totoptdep(*,0,0))-1 do begin - for j=0, n_elements(totoptdep(0,*,0))-1 do begin - for each=0,nday-1 do begin - subs=each*12 + which_utc_is_TES - ref(each)=totoptdep(i,j,subs) - subsarr(each)=subs - endfor - ;;yeah=interpol(ref,subsarr,findgen(ntime), /spline) - yeah=interpol(ref,subsarr,findgen(ntime) ) - totoptdep(i,j,*)=yeah - endfor -endfor - -;;;********************************************************************* - - - - - -;-------------------------------- -print,'writing netcdf file' - -filencdf='input_totoptdep_tot.nc' -cdfid=ncdf_create(filencdf,/clobber) - -londimid=ncdf_dimdef(cdfid,'lon',lon) -latdimid=ncdf_dimdef(cdfid,'lat',lat) -tdimid=ncdf_dimdef(cdfid,'time',ntime) -lonid=ncdf_vardef(cdfid,'lon',[londimid],/float) -ncdf_attput,cdfid,lonid,'long_name','longitude' -ncdf_attput,cdfid,lonid,'units','degrees_east' -latid=ncdf_vardef(cdfid,'lat',[latdimid],/float) -ncdf_attput,cdfid,latid,'long_name','latitude' -ncdf_attput,cdfid,latid,'units','degrees_north' -timeid=ncdf_vardef(cdfid,'time',[tdimid],/float) -ncdf_attput,cdfid,timeid,'long_name','model time' -ncdf_attput,cdfid,timeid,'units','hours since 00:00:00' -varid1=ncdf_vardef(cdfid,fieldname1,[londimid,latdimid,tdimid],/float) -ncdf_attput,cdfid,varid1,'Physics_diagnostic',outfieldstring1 -ncdf_attput,cdfid,varid1,'units',unit1 -ncdf_control,cdfid,/endef - -ncdf_varput,cdfid,lonid,longi -ncdf_varput,cdfid,latid,lati -ncdf_varput,cdfid,timeid,timear -ncdf_varput,cdfid,varid1,totoptdep - -ncdf_close,cdfid - -end Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/fix_no_info_OXF.inc =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/fix_no_info_OXF.inc (revision 1439) +++ (revision ) @@ -1,16 +1,0 @@ - - -!!OXFORD FIX -PRINT *, "prescribed day ..." -!param(4)=1081.-669. -!param(4)=2041.-669.*3. -!param(4)=1021.-669. -!param(4)=1051.-669. -!param(4)=1021. -param(4)=1051. -!param(4)=1201. -!!! -!!! -interval=2 -!!OXFORD FIX - Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/mountain_new.nc_info =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/mountain_new.nc_info (revision 1439) +++ (revision ) @@ -1,18 +1,0 @@ -netcdf mountain_new { -dimensions: - lon = 72 ; - lat = 36 ; - time = 1 ; -variables: - float lon(lon) ; - lon:long_name = "longitude" ; - lon:units = "degrees_east" ; - float lat(lat) ; - lat:long_name = "latitude" ; - lat:units = "degrees_north" ; - float time(time) ; - time:long_name = "model time" ; - time:units = "hours since 00:00:00" ; - float orography(lat, lon) ; - orography:long_name = "orography" ; -} Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/smoothoxf.pro =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/smoothoxf.pro (revision 1439) +++ (revision ) @@ -1,73 +1,0 @@ -pro smoothoxf - -;------------------------------------ -file='input_totoptdep.nc' -smoothampl=4 -;------------------------------------ - -cdfid = ncdf_open(file) -; -; coordinates -; -varid=ncdf_varid(cdfid,'lat') & ncdf_varget, cdfid, varid, lati & lat=n_elements(lati) -varid=ncdf_varid(cdfid,'lon') & ncdf_varget, cdfid, varid, longi & lon=n_elements(longi) -varid=ncdf_varid(cdfid,'time') & ncdf_varget, cdfid, varid, time & ntime=n_elements(time) -; -; field -; -varid=ncdf_varid(cdfid,'totoptdep') & ncdf_varget, cdfid, varid, totoptdep -help, totoptdep -totoptdep=smooth(totoptdep,[smoothampl,smoothampl,0],/EDGE_TRUNCATE) ;; SMOOTH wrt SPACE -fieldname='totoptdep' -outfieldstring='Total dust optical depth (tauref/700Pa)' -unit='Opacity/Pa' - -goto, nointerp -;;;********************************************************************* -; -; INTERPOLATE wrt TIME -; -which_utc_is_TES=3 -;; for Hellas 8AM is 02PM LT, i.e. TES measurements -;; IDL convention, starts 0 !!!! -nday=floor(float(ntime-1)/12.)-1 & ref=fltarr(nday) & subsarr=fltarr(nday) -for i=0, n_elements(totoptdep(*,0,0))-1 do begin - for j=0, n_elements(totoptdep(0,*,0))-1 do begin - for each=0,nday-1 do begin - subs=each*12 + which_utc_is_TES - ref(each)=totoptdep(i,j,subs) - subsarr(each)=subs - endfor - ;;yeah=interpol(ref,subsarr,findgen(ntime), /spline) - yeah=interpol(ref,subsarr,findgen(ntime) ) - totoptdep(i,j,*)=yeah - endfor -endfor -;;;********************************************************************* -nointerp: - - -;-------------------------------- -print,'writing netcdf file' - -filencdf='input_totoptdep_tot.nc' & cdfid=ncdf_create(filencdf,/clobber) -;; -londimid=ncdf_dimdef(cdfid,'lon',lon) & latdimid=ncdf_dimdef(cdfid,'lat',lat) & tdimid=ncdf_dimdef(cdfid,'time',ntime) -;; -lonid=ncdf_vardef(cdfid,'lon',[londimid],/float) & ncdf_attput,cdfid,lonid,'long_name','longitude' & ncdf_attput,cdfid,lonid,'units','degrees_east' -;; -latid=ncdf_vardef(cdfid,'lat',[latdimid],/float) & ncdf_attput,cdfid,latid,'long_name','latitude' & ncdf_attput,cdfid,latid,'units','degrees_north' -;; -timeid=ncdf_vardef(cdfid,'time',[tdimid],/float) & ncdf_attput,cdfid,timeid,'long_name','model time' & ncdf_attput,cdfid,timeid,'units','hours since 00:00:00' -;; -varid=ncdf_vardef(cdfid,fieldname,[londimid,latdimid,tdimid],/float) -ncdf_attput,cdfid,varid,'Physics_diagnostic',outfieldstring -ncdf_attput,cdfid,varid,'units',unit -;; -ncdf_control,cdfid,/endef -;; -ncdf_varput,cdfid,lonid,longi & ncdf_varput,cdfid,latid,lati & ncdf_varput,cdfid,timeid,time & ncdf_varput,cdfid,varid,totoptdep -;; -ncdf_close,cdfid - -end Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/compile_newdiag =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/compile_newdiag (revision 1440) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/compile_newdiag (revision 1440) @@ -0,0 +1,79 @@ +#!/bin/sh + +# ======================================================= +# Compiler script for newdiag.F90 interpol_soil.F90 +# Liam Steele June 2014 +# ======================================================= + +# Determine which Fortran compiler to use +fcomp=0 +echo "Which Fortran compiler do you want to use?" +echo "1. ifort" +echo "2. g95" +echo "3. pgf90" +while [ $fcomp != 1 -a $fcomp != 2 -a $fcomp != 3 ] ; do + read fcomp + if [ $fcomp != 1 -a $fcomp != 2 -a $fcomp != 3 ] ; then + echo "Not a valid compiler choice. Please try again:" + fi +done + +# Look for NetCDF libraries +if [ -n "$NETCDF" ] ; then + echo "Will use NetCDF in directory: $NETCDF" +else + success="" + while [ -z "$success" ] ; do + echo "Cannot find NetCDF library. Please enter full path:" + read ncdfpath + if [ ! -d "$ncdfpath" ] ; then + echo "invalid path: $ncdfpath" ; continue + success="" + else + success="yes" + fi + netcdfipath=$ncdfpath + done + echo "For future use you should put this path into your bash profile, i.e." + echo "export NETCDF="$ncdfpath +fi + +# Check diagfi.nc exists for conversion +if ! [ -e "diagfi.nc" ] ; then + echo 'Please link a valid diagfi.nc file' + exit +fi + +# Check we are in the ukv_newdiag directory +thispwd=`echo $(pwd)` +length=${#thispwd} +substr=`echo $thispwd | cut -c$((length-10))-$((length))` +if [ $substr != "ukv_newdiag" ] ; then + echo 'Not in ukv_newdiag directory' + exit +fi + +# Since all checks are complete, create new diagfi. If you want to use a +# compiler other than ifort you'll have to change the compiler flags +echo 'Compiling executable...' +if [ -e "newdiag.exe" ] ; then + rm newdiag.exe +fi +if [ $fcomp == 1 ] ; then + ifort -convert big_endian -I$NETCDF newdiag.F90 interpol_soil.F90 -o newdiag.exe -L$NETCDF -lnetcdf +elif [ $fcomp == 2 ] ; then + g95 -fno-second-underscore -I$NETCDF newdiag.F90 interpol_soil.F90 -o newdiag.exe -L$NETCDF -lnetcdf +elif [ $fcomp == 3 ] ; then + pgf90 -byteswapio -I$NETCDF newdiag.F90 interpol_soil.F90 -o newdiag.exe -L$NETCDF -lnetcdf +fi +chmod u+x newdiag.exe +./newdiag.exe + +# Create link in PREP_MARS +echo 'Linking new diagfi to PREP_MARS' +cd .. +ln -sf ukv_newdiag/diagfi_new.nc input_diagfi.nc + +echo 'Done!' + +exit Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/interpol_soil.F90 =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/interpol_soil.F90 (revision 1440) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/interpol_soil.F90 (revision 1440) @@ -0,0 +1,114 @@ +! ==================================================================== +! +! This subroutine converts arrays written on regolith levels to arrays +! on soil levels (sometimes different layers had to be specified for +! the regolith scheme (compared to the soil layers) to ensure numerical +! stability at high obliquities) +! +! Liam Steele March 2015 +! +! ==================================================================== + +subroutine interpol_soil(mlayer,mlayer_r,nsoil,nregmx,id,old_array,new_array) +implicit none + +integer :: ik, jk, it, ii, ij +integer :: nsoil, nregmx +real :: lay1, alpha +integer, dimension(4) :: id +real, dimension(nsoil) :: mlayer(0:nsoil-1), layer(1:nsoil) +real, dimension(nreg) :: mlayer_r(0:nregmx-1), layer_r(1:nregmx) +real, dimension(id(1),id(2),id(3),id(4)) :: old_array +real, dimension(id(1),id(2),nsoil,id(4)) :: new_array +integer rloop +real :: lb, ub, lb_r, ub_r +real :: frac(nsoil,nregmx) +integer :: frac_layers(nsoil,nregmx) +integer :: frac_count(nsoil) +integer :: bl(nregmx) + +! Calculate soil lower layer boundaries +lay1 = sqrt(mlayer(0)*mlayer(1)) +alpha = mlayer(1)/mlayer(0) +do jk=1,nsoil + layer(jk) = lay1*(alpha**(jk-1)) +enddo + +! Calculate regolith lower layer boundaries +lay1 = 0.5e-2 +if (nregmx.eq.20) then + alpha = 2.2 +else if (nregmx.eq.30) then + alpha = 1.44 +else if (nregmx.eq.40) then + alpha = 1.26 +else if (nregmx.eq.50) then + alpha = 1.18 +else + print*, 'nregmx is not equal to a specified value' + stop +endif +layer_r(1) = lay1 +do jk=2,10 + layer_r(jk) = layer_r(jk-1) + lay1 +enddo +do jk=11,nregmx + layer_r(jk) = layer_r(jk-1) + alpha*(layer_r(jk-1)-layer_r(jk-2)) +enddo + +! Find which soil layer each regolith layer lies within +bl(:) = 0 +do ik=0,nregmx-1 + do jk=0,nsoil-1 + if (mlayer_r(ik).ge.mlayer(jk)) bl(ik+1) = jk+1 + enddo +enddo + +! Calculate which regolith layers overlap each soil layer, and their contributing fraction +do jk=1,nsoil + lb = layer(jk) + if (jk.eq.1) ub = 0.0 + if (jk.gt.1) ub = layer(jk-1) + frac_count(jk) = 0 + + do ik=1,nregmx + lb_r = layer_r(ik) + if (ik.eq.1) ub_r = 0.0 + if (ik.gt.1) ub_r = layer_r(ik-1) + + if (ub.ge.ub_r .and. lb.le.lb_r) then + ! Whole soil layer is inside regolith layer + frac_count(jk) = frac_count(jk) + 1 + frac(jk,frac_count(jk)) = 1.0 + frac_layers(jk,frac_count(jk)) = ik + else if (ub_r.ge.ub .and. lb_r.le.lb) then + ! Whole regolith layer is inside soil layer + frac_count(jk) = frac_count(jk) + 1 + frac(jk,frac_count(jk)) = (lb_r-ub_r)/(lb-ub) + frac_layers(jk,frac_count(jk)) = ik + else + if (lb_r.gt.ub .and. ub_r.le.ub) then + ! Upper overlap + frac_count(jk) = frac_count(jk) + 1 + frac(jk,frac_count(jk)) = (lb_r-ub)/(lb-ub) + frac_layers(jk,frac_count(jk)) = ik + else if (lb_r.gt.lb .and. ub_r.le.lb) then + ! Lower overlap + frac_count(jk) = frac_count(jk) + 1 + frac(jk,frac_count(jk)) = (lb-ub_r)/(lb-ub) + frac_layers(jk,frac_count(jk)) = ik + endif + endif + enddo +enddo + +! Convert array from regolith layers to soil layers +new_array(:,:,:,:) = 0.0 +do ik=0,nsoil-1 + do rloop=1,frac_count(ik+1) + new_array(:,:,ik+1,:) = new_array(:,:,ik+1,:) + old_array(:,:,frac_layers(ik+1,rloop),:)*frac(ik+1,rloop)/sum(frac(ik+1,1:frac_count(ik+1))) + enddo +enddo + +return +end Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/newdiag.F90 =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/newdiag.F90 (revision 1440) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/ukdiag/newdiag.F90 (revision 1440) @@ -0,0 +1,405 @@ +! ==================================================================== +! +! This routine reads in a diagfi.nc file output from the UK-LMD MGCM +! and converts it into a format suitable for producing initial and +! boundary conditions for the mesoscale model. Specifically, the +! arrays are interpolated onto latitudes ranging from +90 --> -90 +! and longitudes ranging from -180 --> +180. The surface geopotential +! is also added if it doesn't already exist, and a check is made to +! see if the starting sol is present in the controle array. +! +! Arrays written on regolith layers are now interpolated back on to +! soil layers for use in the mesoscale model (sometimes different +! layers had to be specified for the regolith scheme (compared to +! the soil layers) to ensure numerical stability at high obliquities). +! +! The new file created is called diagfi_new.nc, which can be linked +! to the PREP_MARS directory as input_diagfi.nc +! +! Liam Steele June 2014, March 2015 +! +! ==================================================================== + +program newdiag +implicit none + +include "netcdf.inc" + +integer :: i,j,k,ierr,nid,nid2,nid3 +integer :: ndims,nvars,ngatts,unlimdimid +integer :: lonid,latid,sigid,soilid,soildim,regid,timeid,contid,phisid +integer :: nlatnew,nlonnew,nsoil,nreg,npos +integer :: vtype,vatts,vdims,solstart +integer, dimension(4) :: dimids +integer, dimension(:), allocatable :: id,dimid,dimarr,loc,corners,edges +real, dimension(:), allocatable :: lat,lon,sig,soil,regol,time,controle +real, dimension(:), allocatable :: newlat,newlon,weight +real, dimension(:,:), allocatable :: geopot,geopot_new,array2,array2_new +real, dimension(:,:,:), allocatable :: array3,array3_new +real, dimension(:,:,:,:), allocatable :: array4,array4_new +real, dimension(:,:,:,:), allocatable :: new_array +character*20, dimension(:), allocatable :: dimname +character*20 :: vname +character*2 :: printid,res +logical :: wphis,skipv + +! Open file for reading +ierr = nf_open("diagfi.nc",nf_nowrite,nid) +if (ierr /= nf_noerr) then + print*, 'Oh dear: cannot find input_diagfi.nc' + stop +endif + +! Check properties of diagfi +ierr = nf_inq(nid, ndims, nvars, ngatts, unlimdimid) +print*, '--------------' +print*, 'File contains:' +write(printid,'(i2)'),ndims +print*, printid//' dimensions' +write(printid,'(i2)'),nvars +print*, printid//' variables' +print*, '--------------' +allocate(dimarr(ndims)) +allocate(dimname(ndims)) + +! Get dimensions +do i = 1, ndims + ierr = nf_inq_dim(nid,i,dimname(i),dimarr(i)) + if (ierr /= nf_noerr) then + print*, 'Whoops: problem reading dimensions in diagfi.nc' + stop + endif + if (trim(dimname(i)) == 'latitude' .or. trim(dimname(i)) == 'lat') then + dimname(i) = 'latitude' + allocate(lat(dimarr(i))) + ierr = nf_get_var_real(nid,i,lat) + nlatnew = dimarr(i)+1 + else if (trim(dimname(i)) == 'longitude' .or. trim(dimname(i)) == 'lon') then + dimname(i) = 'longitude' + allocate(lon(dimarr(i))) + ierr = nf_get_var_real(nid,i,lon) + nlonnew = dimarr(i)+1 + else if (trim(dimname(i)) == 'sigma') then + allocate(sig(dimarr(i))) + ierr = nf_get_var_real(nid,i,sig) + else if (trim(dimname(i)) == 'soildepth' .or. trim(dimname(i)) == 'soil') then + dimname(i) = 'soildepth' + allocate(soil(dimarr(i))) + ierr = nf_get_var_real(nid,i,soil) + nsoil = dimarr(i) + soildim = i + else if (trim(dimname(i)) == 'regdepth') then + dimname(i) = 'regdepth' + allocate(regol(dimarr(i))) + ierr = nf_get_var_real(nid,i,regol) + nreg = dimarr(i) + else if (trim(dimname(i)) == 'Time' .or. trim(dimname(i)) == 'time') then + dimname(i) = 'Time' + allocate(time(dimarr(i))) + ierr = nf_get_var_real(nid,i,time) + else if (trim(dimname(i)) == 'lentable') then + allocate(controle(dimarr(i))) + ierr = nf_get_var_real(nid,i,controle) + endif +enddo + +! Check for starting sol in diagfi +if (controle(4) == 0) then + print*, 'No starting sol sumber in contole array. Please enter [1-669]:' + read*, solstart + do while (solstart > 669 .or. solstart < 0) + if (solstart > 669) print*, 'Sol > 669. Have another go...' + if (solstart < 0) print*, 'A negative sol is just silly. Try again...' + read*, solstart + enddo + controle(4) = solstart +else if (controle(4) > 669) then + controle(4) = mod(controle(4),669.) +endif + +! Add soil mid-layer depths if not already present +if (soil(1) == 0 .or. soil(1) == 1) then + print*, 'Adding the following soil levels (depths in m):' + do i = 1, nsoil + soil(i) = 2.e-4*(2.**(i-1.5)) + print*, i, soil(i) + enddo +endif + +! Make new lat/lon arrays +allocate(newlat(nlatnew)) +allocate(newlon(nlonnew)) +do i = 1, nlatnew + newlat(i) = 90. + (i-1)*(lat(2)-lat(1)) +enddo +do i = 1, nlonnew + newlon(i) = -180. + (i-1)*(lon(2)-lon(1)) +enddo + +! Set weights and locations for latitudinal interpolation (assumes linear spacing between lats) +allocate(weight(nlatnew)) +allocate(loc(nlatnew)) +weight(:) = 0.5 +weight(1) = 1.5 +weight(nlatnew) = -0.5 +do i = 1, nlatnew + loc(i) = i-1 +enddo +loc(1) = 1 +loc(nlatnew) = nlatnew-2 + +! Create new netcdf file to write data to (doesn't matter about writing descriptions and units) +ierr = nf_create('diagfi_new.nc',nf_clobber,nid2) + +! Define dimensions +do i = 1, ndims + write(printid,'(i2)'),i + print*, printid//' writing '//trim(dimname(i))//' to file' + if (trim(dimname(i)) == 'longitude') then + ierr = nf_def_dim(nid2,'longitude',nlonnew,lonid) + ierr = nf_def_var(nid2,'longitude',nf_float,1,lonid,lonid) + else if (trim(dimname(i)) == 'latitude') then + ierr = nf_def_dim(nid2,'latitude',nlatnew,latid) + ierr = nf_def_var(nid2,'latitude',nf_float,1,latid,latid) + else if (trim(dimname(i)) == 'sigma') then + ierr = nf_def_dim(nid2,'sigma',size(sig),sigid) + ierr = nf_def_var(nid2,'sigma',nf_float,1,sigid,sigid) + else if (trim(dimname(i)) == 'soildepth') then + ierr = nf_def_dim(nid2,'soildepth',size(soil),soilid) + ierr = nf_def_var(nid2,'soildepth',nf_float,1,soilid,soilid) + else if (trim(dimname(i)) == 'regdepth') then + ierr = nf_def_dim(nid2,'regdepth',size(regol),regid) + ierr = nf_def_var(nid2,'regdepth',nf_float,1,regid,regid) + else if (trim(dimname(i)) == 'Time') then + ierr = nf_def_dim(nid2,'Time',size(time),timeid) + ierr = nf_def_var(nid2,'Time',nf_float,1,timeid,timeid) + else if (trim(dimname(i)) == 'lentable') then + ierr = nf_def_dim(nid2,'lentable',size(controle),contid) + ierr = nf_def_var(nid2,'controle',nf_float,1,contid,contid) + endif + if (ierr /= nf_noerr) then + print*, "Yikes! Problem defining dimensions in new_diagfi.nc" + stop + endif +enddo + +ierr = nf_enddef(nid2) + +! Write dimension arrays +ierr = nf_put_var_real(nid2,lonid,newlon) +ierr = nf_put_var_real(nid2,latid,newlat) +ierr = nf_put_var_real(nid2,sigid,sig) +ierr = nf_put_var_real(nid2,soilid,soil) +ierr = nf_put_var_real(nid2,regid,regol) +ierr = nf_put_var_real(nid2,timeid,time) +ierr = nf_put_var_real(nid2,contid,controle) + +! Read each variable in turn from original file, convert to new lat/lon and write to new file +npos = ndims+1 +wphis = .true. +do i = ndims+1, nvars + + ! Read variable dimensions etc. + ierr = nf_inq_var(nid, i, vname, vtype, vdims, dimids, vatts) + write(printid,'(i2)'),i + if (trim(vname) == 'phisinit') wphis = .false. + allocate(id(vdims)) + allocate(corners(vdims)) + allocate(edges(vdims)) + allocate(dimid(vdims)) + + ! Require lon/lat to be first two dimensions (everything should satisfy this condition, but you never know!) + skipv = .true. + if (vdims >= 2) then + if (dimname(dimids(1)) == 'longitude' .and. dimname(dimids(2)) == 'latitude') skipv = .false. + if (dimname(dimids(1)) == 'latitude' .and. dimname(dimids(2)) == 'longitude') skipv = .false. + if (skipv) then + print*, " > "//trim(vname)//": lat and lon aren't first two dimensions. Skipping" + cycle + endif + endif + + ! Change certain names to match LMD (more might need added later) + if (trim(vname) == 'h2ommr') vname = 'q02' + if (trim(vname) == 'icemmr') vname = 'q01' + if (trim(vname) == 'h2o_ice_s') vname = 'qsurf01' + + print*, printid//' writing '//trim(vname)//' to file' + + ! Define arrays for specifying location in diagfi + do k = 1, vdims + id(k) = dimarr(dimids(k)) + dimid(k) = dimids(k) + corners(k) = 1 + edges(k) = dimarr(dimids(k)) + enddo + edges(lonid) = nlonnew + edges(latid) = nlatnew + edges(vdims) = 1 + + ! Allocate and read old arrays + if (vdims == 2) then + allocate(array2(id(1),id(2))) + ierr = nf_get_var_real(nid,i,array2) + else if (vdims == 3) then + allocate(array3(id(1),id(2),id(3))) + ierr = nf_get_var_real(nid,i,array3) + else if (vdims == 4) then + allocate(array4(id(1),id(2),id(3),id(4))) + ierr = nf_get_var_real(nid,i,array4) + endif + + ! Check if regolith depth is the vertical dimension (if so, will convert fields to soil depths) + do k = 1, vdims + if (dimname(dimids(k)).eq.'regdepth') then + + if (vdims == 3) then + allocate(new_array(id(1),id(2),dimarr(soildim),1)) + call interpol_soil(soil,regol,nsoil,nreg,(/id(1:vdims),1/),array3,new_array) + deallocate(array3) + allocate(array3(id(1),id(2),dimarr(soildim))) + array3 = new_array(1:id(1),1:id(2),1:dimarr(soildim),1) + deallocate(new_array) + else if (vdims == 4) then + allocate(new_array(id(1),id(2),dimarr(soildim),id(4))) + call interpol_soil(soil,regol,nsoil,nreg,id(1:vdims),array4,new_array) + deallocate(array4) + allocate(array4(id(1),id(2),dimarr(soildim),id(4))) + array4 = new_array + deallocate(new_array) + endif + + dimids(k) = soildim + dimid(k) = soildim + id(k) = dimarr(dimids(k)) + edges(k) = dimarr(dimids(k)) + + exit + endif + enddo + + ! Allocate new arrays + id(lonid) = nlonnew + id(latid) = nlatnew + if (vdims == 2) allocate(array2_new(id(1),id(2))) + if (vdims == 3) allocate(array3_new(id(1),id(2),id(3))) + if (vdims == 4) allocate(array4_new(id(1),id(2),id(3),id(4))) + + ! Transform onto new lat-lon grid (loops assume array dimensions ordered lon, lat, sig, time) + if (vdims == 2) then + do k = 1, nlatnew + array2_new(1:nlonnew-1,k) = weight(k)*array2(:,loc(k)) + (1-weight(k))*array2(:,loc(k)+1) + enddo + array2_new(nlonnew,:) = array2_new(1,:) + array2_new(:,1) = sum(array2_new(:,2))/nlonnew + array2_new(:,nlatnew) = sum(array2_new(:,nlatnew-1))/nlonnew + if (minval(array2) > 0.0) then + where (array2_new < 0.0) array2_new = 0.0 + endif + else if (vdims == 3) then + do k = 1, nlatnew + array3_new(1:nlonnew-1,k,:) = weight(k)*array3(:,loc(k),:) + (1-weight(k))*array3(:,loc(k)+1,:) + enddo + array3_new(nlonnew,:,:) = array3_new(1,:,:) + do k = 1, id(3) + array3_new(:,1,k) = sum(array3_new(:,2,k))/nlonnew + array3_new(:,nlatnew,k) = sum(array3_new(:,nlatnew-1,k))/nlonnew + enddo + if (minval(array3) > 0.0) then + where (array3_new < 0.0) array3_new = 0.0 + endif + else if (vdims == 4) then + do k = 1, nlatnew + array4_new(1:nlonnew-1,k,:,:) = weight(k)*array4(:,loc(k),:,:) + (1-weight(k))*array4(:,loc(k)+1,:,:) + enddo + array4_new(nlonnew,:,:,:) = array4_new(1,:,:,:) + do k = 1, id(4) + do j = 1, id(3) + array4_new(:,1,j,k) = sum(array4_new(:,2,j,k))/nlonnew + array4_new(:,nlatnew,j,k) = sum(array4_new(:,nlatnew-1,j,k))/nlonnew + enddo + enddo + if (minval(array4) > 0.0) then + where (array4_new < 0.0) array4_new = 0.0 + endif + else + print*, ' > '//trim(vname)//': code does not interpolate 1D variables' + endif + + ! Write variable to new file + ierr = nf_redef(nid2) + ierr = nf_def_var(nid2,trim(vname),nf_float,vdims,dimid,npos) + + ierr = nf_enddef(nid2) + if (vdims == 2) then + ierr = nf_put_var_real(nid2,npos,array2_new) + else + do k = 1, id(vdims) + corners(vdims) = k + if (vdims == 3) ierr = nf_put_vara_real(nid2,npos,corners,edges,array3_new(:,:,k)) + if (vdims == 4) ierr = nf_put_vara_real(nid2,npos,corners,edges,array4_new(:,:,:,k)) + enddo + if (vdims == 3) deallocate(array3,array3_new) + if (vdims == 4) deallocate(array4,array4_new) + endif + + if (ierr /= nf_noerr) then + print*, "Oh no! Problem writing "//trim(vname) + stop + endif + + deallocate(id,corners,edges,dimid) + npos = npos + 1 + +enddo + +! Add geopotential height to file +if (wphis) then + print*, 'Adding phisinit to file' + allocate(id(2)) + + if (nlatnew-1 == 8) res = '5' + if (nlatnew-1 == 14) res = '10' + if (nlatnew-1 == 24) res = '21' + if (nlatnew-1 == 36) res = '31' + if (nlatnew-1 == 48) res = '42' + if (nlatnew-1 == 72) res = '63' + if (nlatnew-1 == 96) res = '85' + if (nlatnew-1 == 144) res = '127' + if (nlatnew-1 == 192) res = '170' + allocate(geopot(nlonnew-1,nlatnew-1)) + allocate(geopot_new(nlonnew,nlatnew)) + ierr = nf_open("phisinit_files/phisinit-T"//trim(res)//".nc",nf_nowrite,nid3) + ierr = nf_inq_varid(nid3,"phis",phisid) + if (ierr /= nf_noerr) then + print*, "Oh no! Can't find geopotential file. Please link from the ukv_phisinit folder:" + print*, "ln -sf $MMM/your_comp_dir/PREP_MARS/ukv_phisinit/phisinit-T"//trim(res)//".nc ." + stop + endif + ierr = nf_get_var_real(nid3,phisid,geopot) + ierr = nf_close(nid3) + + do k = 1, nlatnew + geopot_new(1:nlonnew-1,k) = weight(k)*geopot(:,loc(k)) + (1-weight(k))*geopot(:,loc(k)+1) + enddo + geopot_new(nlonnew,:) = geopot_new(1,:) + geopot_new(:,1) = sum(geopot_new(:,2))/nlonnew + geopot_new(:,nlatnew) = sum(geopot_new(:,nlatnew-1))/nlonnew + + id(1) = lonid + id(2) = latid + ierr = nf_redef(nid2) + ierr = nf_def_var(nid2,"phisinit",nf_float,2,id,npos) + ierr = nf_enddef(nid2) + ierr = nf_put_var_real(nid2,npos,geopot_new) + if (ierr /= nf_noerr) then + print*, "Oh no! Problem writing phisinit" + stop + endif +endif + +print*, 'Closing files...' +ierr = nf_close(nid) +ierr = nf_close(nid2) + +end Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/makegcm_mpifort =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/makegcm_mpifort (revision 1439) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/makegcm_mpifort (revision 1440) @@ -183,9 +183,9 @@ #-ffixed-form # -set optim=" -O2 -static -zero -align commons" +set optim=" -O2 -static -zero -align commons -mcmodel=large -shared-intel" #-ffree-line-length-huge" -set optim90=" -O2 -static -zero -align commons" +set optim90=" -O2 -static -zero -align commons -mcmodel=large -shared-intel" #-ffree-line-length-huge" -set optimtru90=" -O2 -static -zero -align commons" +set optimtru90=" -O2 -static -zero -align commons -mcmodel=large -shared-intel" #-ffree-line-length-huge" set opt_link=" -L$NCDFLIB -lnetcdf" Index: trunk/MESOSCALE/LMD_MM_MARS/makemeso =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/makemeso (revision 1439) +++ trunk/MESOSCALE/LMD_MM_MARS/makemeso (revision 1440) @@ -36,4 +36,5 @@ scenario="" from_scratch=0 + while getopts "drc:njhgpfs:xe" options; do case $options in @@ -492,4 +493,7 @@ sed s+"-cc=icc"+" "+g configure.wrf > yeah ; mv -f yeah configure.wrf sed s+"-DIFORT_KLUDGE"+" "+g configure.wrf > yeah ; mv -f yeah configure.wrf + sed s+"-O1"+"-O1 -mcmodel=medium -shared-intel"+g configure.wrf > yeah ; mv -f yeah configure.wrf + sed s+"-O2"+"-O2 -mcmodel=medium -shared-intel"+g configure.wrf > yeah ; mv -f yeah configure.wrf + sed s+"-O3"+"-O3 -mcmodel=medium -shared-intel"+g configure.wrf > yeah ; mv -f yeah configure.wrf if [ ${debug} -ne 0 ] # not working for xlf! then