Changeset 2401
- Timestamp:
- Jul 7, 2020, 3:49:34 PM (4 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r2400 r2401 3124 3124 == 07/07/2020 == DB 3125 3125 Add kstar parameter to control kmin value. kmi=1/lambda_max, to ensure the "subgrid scale" characteristic, we have to constrain the GW's wavelength by the size of the mesh 3126 3127 == 07/07/2020 == EM 3128 Improve "zrecast" utility: add possibility to transfer 3D (lon-lat-time) 3129 variables to output file and add some extra sanity checks of user inputs. -
trunk/LMDZ.MARS/util/zrecast.F90
r2140 r2401 79 79 integer nbvarinfile ! # of variables in input file 80 80 integer nbattr ! # of attributes of a given variable in input file 81 integer nbvar3dinfile ! # of 3D (lon,lat,time) variables in input file 81 82 integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file 82 83 integer outfid ! NetCDF output file ID … … 97 98 character (len=64), dimension(:), allocatable :: var 98 99 ! var(): names of variables that will be processed 100 integer, allocatable :: varndims(:) ! number of dimensions of var(:) variables 99 101 integer nbvar ! # of variables to process 100 102 integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile) … … 120 122 real,dimension(:,:,:,:),allocatable :: ra_gcm ! radial distance (m) to 121 123 ! center of Mars, (of GCM grid points) 124 real,dimension(:,:,:),allocatable :: surfdata ! to store a GCM 3D dataset 122 125 real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset 123 126 real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset … … 186 189 write(*,*)" The following variables have been found:" 187 190 nbvar4dinfile=0 191 nbvar3dinfile=0 188 192 do i=1,nbvarinfile 189 193 ! get name of variable # i … … 195 199 write(*,*) trim(tmpvarname) 196 200 endif 201 if (tmpndims.eq.3) then 202 nbvar3dinfile=nbvar3dinfile+1 203 write(*,*) trim(tmpvarname) 204 endif 197 205 enddo 198 206 199 allocate(var(nbvar4dinfile ),stat=ierr)207 allocate(var(nbvar4dinfile+nbvar3dinfile),stat=ierr) 200 208 if (ierr.ne.0) then 201 write(*,*) "Error! failed memory allocation of var(nbvar4dinfile)" 209 write(*,*) "Error! failed memory allocation of var(nbvar4dinfile+nbvar3dinfile)" 210 stop 211 endif 212 allocate(varndims(nbvar4dinfile+nbvar3dinfile),stat=ierr) 213 if (ierr.ne.0) then 214 write(*,*) "Error! failed memory allocation of varndims(nbvar4dinfile+nbvar3dinfile)" 202 215 stop 203 216 endif … … 213 226 ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid) 214 227 if (ierr.eq.NF_NOERR) then ! valid name 215 ! also check that it is indeed a 4D variable228 ! also check that it is indeed a 3D or 4D variable 216 229 ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims) 217 if ( tmpndims.eq.4) then230 if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then 218 231 nbvar=nbvar+1 219 var(nbvar)=tmpvarname 232 var(nbvar)=tmpvarname 233 varndims(nbvar)=tmpndims 220 234 else 221 write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'235 write(*,*) 'Error: ',trim(tmpvarname),' is not a 3D or 4D variable' 222 236 write(*,*) ' (we''ll skip that one)' 223 237 endif … … 233 247 nbvar=0 234 248 do i=1,nbvarinfile 235 ! look for 4D variables249 ! look for 3D and 4D variables 236 250 ierr=NF_INQ_VARNDIMS(infid,i,tmpndims) 237 if ( tmpndims.eq.4) then251 if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then 238 252 nbvar=nbvar+1 239 253 ! get the corresponding name 240 254 ierr=NF_INQ_VARNAME(infid,i,tmpvarname) 241 255 var(nbvar)=tmpvarname 256 varndims(nbvar)=tmpndims 242 257 endif 243 258 enddo … … 622 637 if (.not.auto_mcd_levels) then 623 638 write(*,*) "" 624 write(*,*) "Enter min and max of vertical coordinate (Pa or m)"625 write(*,*) " (in that order and on the same line)"626 639 if (ztype.eq.1) then ! pressure coordinate 627 read(*,*) pmin,pmax 640 write(*,*) "Enter min and max of vertical coordinate (Pa)" 641 write(*,*) " (in that order and on the same line)" 642 read(*,*,iostat=ierr) pmin,pmax 643 if ((ierr/=0).or.(pmin.gt.pmax)) then 644 write(*,*) "Error: wrong vertical coordinate boundaries" 645 write(*,*) " pmin=",pmin 646 write(*,*) " pmax=",pmax 647 stop 648 endif 628 649 else ! altitude coordinate, except in MCD case 629 read(*,*) zamin,zamax 650 write(*,*) "Enter min and max of vertical coordinate (m)" 651 write(*,*) " (in that order and on the same line)" 652 read(*,*,iostat=ierr) zamin,zamax 653 if ((ierr/=0).or.(zamin.gt.zamax)) then 654 write(*,*) "Error: wrong vertical coordinate boundaries" 655 write(*,*) " zamin=",zamin 656 write(*,*) " zamax=",zamax 657 stop 658 endif 630 659 endif 631 660 endif … … 634 663 if (ztype.eq.1) then ! pressure coordinate 635 664 write(*,*) "Number of levels along vertical coordinate?" 636 read(*,*) nblev 665 read(*,*,iostat=ierr) nblev 666 if (ierr/=0) then 667 write(*,*) "Error: bad input value for nblev ", nblev 668 stop 669 endif 637 670 allocate(plevel(nblev),stat=ierr) 638 671 if (ierr.ne.0) then … … 652 685 else if (ztype.eq.2) then ! above areoid heights 653 686 write(*,*) "Number of levels along vertical coordinate?" 654 read(*,*) nblev 687 read(*,*,iostat=ierr) nblev 688 if (ierr/=0) then 689 write(*,*) "Error: bad input value for nblev ", nblev 690 stop 691 endif 655 692 allocate(zareoid(nblev),stat=ierr) 656 693 if (ierr.ne.0) then … … 680 717 else 681 718 write(*,*) "Number of levels along vertical coordinate?" 682 read(*,*) nblev 719 read(*,*,iostat=ierr) nblev 720 if (ierr/=0) then 721 write(*,*) "Error: bad input value for nblev ", nblev 722 stop 723 endif 683 724 allocate(zsurface(nblev),stat=ierr) 684 725 if (ierr.ne.0) then … … 697 738 else ! distance to center of planet 698 739 write(*,*) "Number of levels along vertical coordinate?" 699 read(*,*) nblev 740 read(*,*,iostat=ierr) nblev 741 if (ierr/=0) then 742 write(*,*) "Error: bad input value for nblev ", nblev 743 stop 744 endif 700 745 allocate(zradius(nblev),stat=ierr) 701 746 if (ierr.ne.0) then … … 716 761 write(*,*) "" 717 762 write(*,*) "Number of levels along vertical coordinate?" 718 read(*,*) nblev 763 read(*,*,iostat=ierr) nblev 764 if (ierr/=0) then 765 write(*,*) "Error: bad input value for nblev ", nblev 766 stop 767 endif 719 768 if (ztype.eq.1) then ! pressure coordinate 720 769 allocate(plevel(nblev),stat=ierr) … … 728 777 write(*,*) " (one value per line)" 729 778 do i=1,nblev 730 read(*,*) plevel(i) 779 read(*,*,iostat=ierr) plevel(i) 780 if (ierr/=0) then 781 write(*,*) "Error: bad input value for plevel(i):",plevel(i) 782 stop 783 endif 731 784 enddo 732 785 else if (ztype.eq.2) then ! above areoid heights … … 740 793 write(*,*) " from min to max (one value per line)" 741 794 do i=1,nblev 742 read(*,*) zareoid(i) 795 read(*,*,iostat=ierr) zareoid(i) 796 if (ierr/=0) then 797 write(*,*) "Error: bad input value for zareoid(i):",zareoid(i) 798 stop 799 endif 743 800 enddo 744 801 else if (ztype.eq.3) then ! above local surface … … 752 809 write(*,*) " from min to max (one value per line)" 753 810 do i=1,nblev 754 read(*,*) zsurface(i) 811 read(*,*,iostat=ierr) zsurface(i) 812 if (ierr/=0) then 813 write(*,*) "Error: bad input value for zsurface(i):",zsurface(i) 814 stop 815 endif 755 816 enddo 756 817 else ! distance to center of planet … … 764 825 write(*,*) " from min to max (one value per line)" 765 826 do i=1,nblev 766 read(*,*) zradius(i) 827 read(*,*,iostat=ierr) zradius(i) 828 if (ierr/=0) then 829 write(*,*) "Error: bad input value for zradius(i):",zradius(i) 830 stop 831 endif 767 832 enddo 768 833 endif … … 1317 1382 endif 1318 1383 1319 ! 3.3.3 Define 4D variables1384 ! 3.3.3 Define 3D and 4D variables 1320 1385 1321 1386 ! add pressure or zareoid … … 1381 1446 stop 1382 1447 endif 1383 do i=1,nbvar 1448 do i=1,nbvar ! loop on 3D+4D variables 1449 if (var(i)=="ps") cycle ! skip surface pressure (already processed) 1450 ! other variables 1384 1451 write(*,*) "" 1385 1452 write(*,*) "Creating ",trim(var(i)) 1386 ! define the variable 1387 ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i)) 1388 if (ierr.ne.NF_NOERR) then 1389 write(*,*) 'Error, could not define variable ',trim(var(i)) 1390 write(*,*) NF_STRERROR(ierr) 1391 stop 1453 if (varndims(i).eq.3) then 1454 ! define the variable 1455 ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,3,surfdatashape,var_id(i)) 1456 if (ierr.ne.NF_NOERR) then 1457 write(*,*) 'Error, could not define variable ',trim(var(i)) 1458 write(*,*) NF_STRERROR(ierr) 1459 stop 1460 endif 1461 elseif (varndims(i).eq.4) then 1462 ! define the variable 1463 ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i)) 1464 if (ierr.ne.NF_NOERR) then 1465 write(*,*) 'Error, could not define variable ',trim(var(i)) 1466 write(*,*) NF_STRERROR(ierr) 1467 stop 1468 endif 1469 else 1470 write(*,*) "Unexpected value for varndims(i)!" 1471 write(*,*) " i=",i," varndims(i)=",varndims(i) 1472 stop 1392 1473 endif 1393 1474 … … 1605 1686 endif 1606 1687 1688 ! allocate surfdata() to store input 3D GCM data 1689 allocate(surfdata(lonlength,latlength,timelength),stat=ierr) 1690 if (ierr.ne.0) then 1691 write(*,*) "Error: Failed to allocate surfdata(lonlength,latlength,timelength)" 1692 write(*,*) " lonlength=",lonlength," latlength=",latlength 1693 write(*,*) " timelength=",timelength 1694 stop 1695 endif 1696 1697 do i=1,nbvar ! loop on 3D+4D variables 1698 if (varndims(i).eq.4) cycle ! skip 4D variables 1699 if (var(i)=="ps") cycle ! skip surface pressure (already processed) 1700 ! load the variable from input file 1701 ierr=NF_INQ_VARID(infid,var(i),tmpvarid) 1702 if (ierr.ne.NF_NOERR) then 1703 write(*,*) "Error: Failed to get ",trim(var(i)),"ID" 1704 write(*,*) NF_STRERROR(ierr) 1705 stop 1706 endif 1707 ierr=NF_GET_VAR_REAL(infid,tmpvarid,surfdata) 1708 if (ierr.ne.NF_NOERR) then 1709 write(*,*) "Error: Failed to load ",trim(var(i)) 1710 write(*,*) NF_STRERROR(ierr) 1711 stop 1712 endif 1713 ! write the variable to output file 1714 ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner(1:3),edges(1:3),surfdata) 1715 if (ierr.ne.NF_NOERR) then 1716 write(*,*) "Error: Could not write ",trim(var(i))," data to output file" 1717 write(*,*) NF_STRERROR(ierr) 1718 stop 1719 endif 1720 enddo ! of do i=1,nbvar 1721 1607 1722 !=============================================================================== 1608 1723 ! 4. Interpolate and write 4D variables … … 1629 1744 ! 4.1 If output is in pressure coordinate 1630 1745 if (ztype.eq.1) then 1631 do i=1,nbvar ! loop on 4D variable to process 1746 do i=1,nbvar ! loop on 3D+4D variable to process 1747 ! Skip 3D variables 1748 if (varndims(i).eq.3) cycle 1632 1749 ! identify and read a dataset 1633 1750 ierr=NF_INQ_VARID(infid,var(i),tmpvarid) … … 1675 1792 ! 4.2 If output is in above areoid altitude 1676 1793 if (ztype.eq.2) then 1677 do i=1,nbvar ! loop on 4D variable to process 1794 do i=1,nbvar ! loop on 3D+4D variable to process 1795 ! Skip 3D variables 1796 if (varndims(i).eq.3) cycle 1678 1797 ! identify and read a dataset 1679 1798 ierr=NF_INQ_VARID(infid,var(i),tmpvarid) … … 1726 1845 ! 4.3 If output is in above local surface altitude 1727 1846 if (ztype.eq.3) then 1728 do i=1,nbvar ! loop on 4D variable to process 1847 do i=1,nbvar ! loop on 3D+4D variable to process 1848 ! Skip 3D variables 1849 if (varndims(i).eq.3) cycle 1729 1850 write(*,*) " " 1730 1851 write(*,*) "Processing "//trim(var(i)) … … 1781 1902 ! 4.4 If output is in radial distance to center of planet 1782 1903 if (ztype.eq.4) then 1783 do i=1,nbvar ! loop on 4D variable to process 1904 do i=1,nbvar ! loop on 3D+4D variable to process 1905 ! Skip 3D variables 1906 if (varndims(i).eq.3) cycle 1784 1907 ! identify and read a dataset 1785 1908 ierr=NF_INQ_VARID(infid,var(i),tmpvarid) … … 1836 1959 write(*,*) NF_STRERROR(ierr) 1837 1960 endif 1961 1962 write(*,*) 1963 write(*,*) "zrecast: All is well that ends well." 1838 1964 1839 1965 end program … … 2846 2972 ! data structure of gravity field 2847 2973 double precision v0,omega,ae,gm 2848 2849 2974 double precision clm(0:ndeg,0:ndeg),slm(0:ndeg,0:ndeg) 2975 common /gmm1/v0,omega,ae,gm,clm,slm 2850 2976 integer lmin,lmax 2851 2977 double precision root(nd2p3) 2852 2978 double precision requator 2853 2979 common /sqr/ lmin,lmax,root,requator 2854 2980 2855 2981 !=============================================================================== … … 2865 2991 dlat=lat(jloop) ! double precision version of lat 2866 2992 call geoid(dlon,dlat,dareoid) 2867 rareoid(iloop,jloop)= dareoid2993 rareoid(iloop,jloop)=real(dareoid) 2868 2994 enddo 2869 2995 enddo … … 2901 3027 ! data structure of gravity field 2902 3028 integer ndeg,ndeg2,nd2p3 2903 2904 2905 2906 2907 2908 3029 parameter (ndeg=90,ndeg2=2*ndeg,nd2p3=ndeg2+3) 3030 double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg) 3031 double precision v0,omega,ae,gm !,clm,slm 3032 double precision plm(ndeg) 3033 common /gmm1/ v0,omega,ae,gm,c,s 3034 integer lmin,lmax 2909 3035 double precision root(nd2p3) 2910 3036 double precision r … … 2916 3042 double precision x,xi,sum 2917 3043 2918 2919 2920 3044 ae= 3396000.d0 3045 gm =42828.37d9 3046 omega=0.70882181d-4 2921 3047 ! this value makes the equatorial mean radius equal to 3396 km. 2922 2923 2924 3048 do k=1,nd2p3 3049 root(k)=sqrt(dble(k)) 3050 enddo 2925 3051 ! initialize coefficients c and s 2926 3052 c(:,:)=0 … … 9488 9614 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 9489 9615 9490 9491 9492 9493 9494 9495 9496 9497 9498 9499 9500 ! centrifugal potential 9501 9502 ! 9503 9616 r = 3396000.d0 ! MOLA potential surface 9617 xi = ae/r 9618 m=0 9619 x=0. 9620 sum=0. 9621 call lgndr(lmax,m,x,plm,root) 9622 do l=2, lmax 9623 sum = sum + c(l,m) *xi**l *plm(l-m+1) 9624 enddo 9625 sum = sum+1. 9626 ! centrifugal potential 9627 v0=(gm*sum/r + 0.5*omega**2 * r**2) 9628 ! write(out,*)'v0=', v0 9629 return 9504 9630 9505 9631 end subroutine areoid_ini … … 9508 9634 9509 9635 9510 9636 subroutine geoid(dlon,dlat,rg) 9511 9637 ! marsgeoid.f 9512 9638 ! calculate radii of surface of constant gravity potential on rotating body … … 9521 9647 9522 9648 double precision pi,d2r 9523 9649 parameter (pi=3.141592653589792D0, d2r=pi/180.d0) 9524 9650 ! data structure of gravity field 9525 9651 integer ndeg,ndeg1,ndeg2,nd2p3 9526 9652 parameter (ndeg=90,ndeg1=ndeg+1,ndeg2=2*ndeg,nd2p3=ndeg2+3) 9527 9653 double precision v0,omega,ae,gm 9528 9529 9654 double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg) 9655 common /gmm1/v0,omega,ae,gm,c,s 9530 9656 integer lmin,lmax 9531 9657 double precision root(nd2p3) 9532 9658 double precision r 9533 9534 9535 9659 common /sqr/ lmin,lmax,root,r 9660 double precision plm(ndeg1) 9661 double precision tol 9536 9662 data tol /0.125d0/ ! tolerence on computed value of rg 9537 9663 … … 9541 9667 double precision diff 9542 9668 ! save r 9543 9544 9545 9546 9547 9548 9549 9550 9551 9552 9553 9554 9555 9556 9557 9558 9559 9560 9561 9562 ! centrifugal potential 9563 9669 rg = r 9670 9671 rlon=dlon*d2r 9672 rlat=dlat*d2r 9673 x = sin(rlat) 9674 cslt= cos(rlat) 9675 9676 do i=1,8 ! usually 3 iterations suffice 9677 9678 xi = ae/r 9679 sum = 0. 9680 do m=0, lmax 9681 call lgndr(lmax,m,x,plm,root) 9682 do l=m, lmax 9683 cslm = c(l,m)*cos(m*rlon) + s(l,m)*sin(m*rlon) 9684 sum = sum + cslm*xi**l *plm(l-m+1) 9685 enddo 9686 enddo 9687 sum=sum+1. 9688 ! centrifugal potential 9689 rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0 9564 9690 ! write(out,*) i,r,rg 9565 9691 diff = r-rg 9566 9692 r = rg 9567 9693 if( abs(diff).lt. tol) goto 400 9568 9569 9570 400 9571 9572 9694 enddo ! do i=1,8 9695 9696 400 continue 9697 9698 end subroutine geoid 9573 9699 9574 9700 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 9608 9734 if (LMAX.gt.M+1) then 9609 9735 DO LL=M+2,LMAX 9610 9736 PLL=SQR(2*LL+1)*SQR(2*LL-1)/(SQR(LL+M)*SQR(LL-M)) & 9611 9737 & * (X*PMMP1-PMM*SQR(LL+M-1)*SQR(LL-M-1)/(SQR(2*LL-1)*SQR(2*LL-3))) 9612 9738 PLM(LL-M+1)=OM*PLL
Note: See TracChangeset
for help on using the changeset viewer.