Changeset 14974 for utils/tools_r4.0-HEAD/NESTING
- Timestamp:
- 2021-06-11T10:35:53+02:00 (3 years ago)
- Location:
- utils/tools_r4.0-HEAD/NESTING/src
- Files:
-
- 2 edited
-
agrif_interpolation.f90 (modified) (7 diffs)
-
io_netcdf.f90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
utils/tools_r4.0-HEAD/NESTING/src/agrif_interpolation.f90
r12243 r14974 504 504 REAL*8, POINTER, DIMENSION(:,:) :: lonChild,latChild,latlon_temp => NULL() 505 505 REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL() 506 REAL*8, POINTER, DIMENSION(:,:,:) :: tabinterp3d,tabvar3d => NULL()506 REAL*8, POINTER, DIMENSION(:,:,:) :: tabinterp3d,tabvar3d,tabvar01,tabvar02,tabvar03 => NULL() 507 507 REAL*8, POINTER, DIMENSION(:) :: timedepth_temp,depth => NULL() 508 508 REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL() … … 516 516 LOGICAL :: Interpolation,conservation,Pacifique,Extrapolation,land_level 517 517 ! 518 INTEGER :: deptht,time,i,status,ncid,t,ii,j,nb,numlon,numlat 518 INTEGER :: deptht,time,i,status,ncid,t,ii,j,nb,numlon,numlat,unlimdimid 519 519 520 520 ! 521 521 TYPE(Coordinates) :: G0,G1 522 ! 523 status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid) 524 IF (status/=nf90_noerr) THEN 525 WRITE(*,*)"unable to open netcdf file : ",TRIM(filename) 526 STOP 527 ENDIF 522 528 ! 523 529 !***************** … … 542 548 CALL Read_Ncdf_dim('x',filename,numlon) 543 549 CALL Read_Ncdf_dim('y',filename,numlat) 544 CALL Read_Ncdf_dim('time_counter',filename,time)545 550 IF ( Dims_Existence( 'deptht' , filename ) ) THEN 546 551 CALL Read_Ncdf_dim('deptht',filename,deptht) … … 556 561 deptht = N 557 562 ENDIF 558 563 IF ( Dims_Existence( 'time_counter' , filename ) ) THEN 564 CALL Read_Ncdf_dim('time_counter',filename,time) 565 ! check that time_counter is the unlimited dim 566 status = nf90_inquire(ncid, unlimiteddimid = unlimdimid) 567 IF ( unlimdimid == -1 ) THEN 568 WRITE(*,*)"time_counter should be the unlimited dimension in this file : ",filename 569 WRITE(*,*)" use nco command: ncks -O --mk_rec_dmn time_counter ",filename," ",filename 570 STOP 571 ENDIF 572 ELSEIF( Dims_Existence( 'time' , filename ) ) THEN 573 CALL Read_Ncdf_dim('time',filename,time) 574 ! check that time is the unlimited dim 575 status = nf90_inquire(ncid, unlimiteddimid = unlimdimid) 576 IF ( unlimdimid == -1 ) THEN 577 WRITE(*,*)"time should be the unlimited dimension in this file : ",filename 578 WRITE(*,*)" use nco command: ncks -O --mk_rec_dmn time ",filename," ",filename 579 STOP 580 ENDIF 581 ELSE 582 time = 0 583 ENDIF 584 559 585 ! 560 586 ! retrieve netcdf variable name … … 666 692 IF ( Dims_Existence( 'depthv' , filename ) ) CALL Write_Ncdf_dim('depthv',Child_file,deptht) 667 693 IF ( Dims_Existence( 'depthw' , filename ) ) CALL Write_Ncdf_dim('depthw',Child_file,deptht) 668 IF ( Dims_Existence( 'z' , filename ) ) CALL Write_Ncdf_dim('z',Child_file,deptht) 669 CALL Write_Ncdf_dim('time_counter',Child_file,0) 694 IF ( Dims_Existence( 'z' , filename ) ) CALL Write_Ncdf_dim('z' ,Child_file,deptht) 695 IF ( Dims_Existence( 'time_counter' , filename ) ) CALL Write_Ncdf_dim('time_counter',Child_file,time) 696 IF ( Dims_Existence( 'time' , filename ) ) CALL Write_Ncdf_dim('time' ,Child_file,time) 670 697 671 698 IF( deptht .NE. 1 .AND. deptht .NE. N ) THEN … … 730 757 Child_file,timedepth_temp,'float') 731 758 CALL Copy_Ncdf_att('time_counter',TRIM(filename),Child_file) 759 DEALLOCATE(timedepth_temp) 760 varname = TRIM(Ncdf_varname(i)) 761 Interpolation = .FALSE. 762 ! 763 !copy time from input file to output file 764 CASE('time') 765 ALLOCATE(timedepth_temp(time)) 766 CALL Read_Ncdf_var('time',filename,timedepth_temp) 767 CALL Write_Ncdf_var('time','time', & 768 Child_file,timedepth_temp,'float') 769 CALL Copy_Ncdf_att('time',TRIM(filename),Child_file) 732 770 DEALLOCATE(timedepth_temp) 733 771 varname = TRIM(Ncdf_varname(i)) … … 807 845 END SELECT 808 846 809 ! //////////////// INTERPOLATION FOR 3D VARIABLES /////////////////////////////////////847 ! //////////////// INTERPOLATION FOR 2D VARIABLES ///////////////////////////////////// 810 848 ! 811 IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3) THEN849 IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 2 ) THEN 812 850 ! 813 851 ALLOCATE(detected_pts(numlon,numlat,N)) 814 852 ALLOCATE(masksrc(numlon,numlat)) 815 853 ! 816 ! ******************************LOOP ON TIME******************************************* 817 !loop on time 818 DO t = 1,time 819 ! 820 IF(extrapolation) THEN 821 WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t 822 ELSE 823 WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t 824 ENDIF 825 ! 826 ALLOCATE(tabvar3d(numlon,numlat,1)) 827 ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 828 ! 829 CALL Read_Ncdf_var(varname,filename,tabvar3d,t) 830 ! 831 ! search points where extrapolation is required 832 ! 833 IF(Extrapolation) THEN 834 WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0. 835 IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1) 836 CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar') 837 ELSE 838 masksrc = .TRUE. 839 ENDIF 840 841 IF (t.EQ.1 ) THEN 842 843 SELECT CASE(TRIM(interp_type)) 844 CASE('bilinear') 845 CALL get_remap_matrix(latParent,latChild, & 846 lonParent,lonChild, & 847 masksrc,matrix,src_add,dst_add) 848 849 CASE('bicubic') 850 CALL get_remap_bicub(latParent,latChild, & 851 lonParent,lonChild, & 852 masksrc,matrix,src_add,dst_add) 853 854 END SELECT 855 ! 856 ENDIF 857 ! 858 SELECT CASE(TRIM(interp_type)) 859 CASE('bilinear') 860 CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, & 861 matrix,src_add,dst_add) 862 CASE('bicubic') 863 CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, & 864 matrix,src_add,dst_add) 865 END SELECT 866 ! 867 IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), & 868 G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 869 ! 870 IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1) 871 ! 872 dimnames(1)='x' 873 dimnames(2)='y' 874 dimnames(3)='time_counter' 875 ! 876 CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),& 877 Child_file,tabinterp3d,t,'float') 878 ! 879 DEALLOCATE(tabinterp3d) 880 DEALLOCATE(tabvar3d) 881 !end loop on time 882 END DO 854 IF(extrapolation) THEN 855 WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname) 856 ELSE 857 WRITE(*,*) 'interpolation ',TRIM(varname) 858 ENDIF 859 ! 860 ALLOCATE(tabvar3d(numlon,numlat,1)) 861 ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 862 ! 863 CALL Read_Ncdf_var(varname,filename,tabvar3d) 864 ! 865 ! search points where extrapolation is required 866 ! 867 IF(Extrapolation) THEN 868 WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0. 869 CALL extrap_detect(G0,G1,detected_pts(:,:,1),1) 870 CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar') 871 ELSE 872 masksrc = .TRUE. 873 ENDIF 874 875 SELECT CASE(TRIM(interp_type)) 876 CASE('bilinear') 877 CALL get_remap_matrix(latParent,latChild, & 878 lonParent,lonChild, & 879 masksrc,matrix,src_add,dst_add) 880 881 CASE('bicubic') 882 CALL get_remap_bicub(latParent,latChild, & 883 lonParent,lonChild, & 884 masksrc,matrix,src_add,dst_add) 885 886 END SELECT 887 ! 888 ! 889 SELECT CASE(TRIM(interp_type)) 890 CASE('bilinear') 891 CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, & 892 matrix,src_add,dst_add) 893 CASE('bicubic') 894 CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, & 895 matrix,src_add,dst_add) 896 END SELECT 897 ! 898 IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), & 899 G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 900 ! 901 IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1) 902 ! 903 dimnames(1)='x' 904 dimnames(2)='y' 905 ! 906 CALL Write_Ncdf_var(TRIM(varname),dimnames(1:2),& 907 Child_file,tabinterp3d(:,:,1),'float') 908 ! 909 DEALLOCATE(tabinterp3d) 910 DEALLOCATE(tabvar3d) 883 911 ! 884 912 DEALLOCATE(detected_pts) 885 913 IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add) 886 914 DEALLOCATE( masksrc) 915 916 CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file) 917 ! 918 919 ! //////////////// INTERPOLATION FOR 3D VARIABLES ///////////////////////////////////// 920 ! 921 ELSEIF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3 ) THEN 922 ! 923 IF( DimUnlimited_Var(TRIM(varname),TRIM(filename)) ) THEN 924 925 ALLOCATE(detected_pts(numlon,numlat,N)) 926 ALLOCATE(masksrc(numlon,numlat)) 927 ! 928 ! ******************************LOOP ON TIME******************************************* 929 !loop on time 930 DO t = 1,time 931 ! 932 IF(extrapolation) THEN 933 WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t 934 ELSE 935 WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t 936 ENDIF 937 ! 938 ALLOCATE(tabvar3d(numlon,numlat,1)) 939 ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 940 ! 941 CALL Read_Ncdf_var(varname,filename,tabvar3d,t) 942 ! 943 ! search points where extrapolation is required 944 ! 945 IF(Extrapolation) THEN 946 WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0. 947 IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1) 948 CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar') 949 ELSE 950 masksrc = .TRUE. 951 ENDIF 952 953 IF (t.EQ.1 ) THEN 954 955 SELECT CASE(TRIM(interp_type)) 956 CASE('bilinear') 957 CALL get_remap_matrix(latParent,latChild, & 958 lonParent,lonChild, & 959 masksrc,matrix,src_add,dst_add) 960 961 CASE('bicubic') 962 CALL get_remap_bicub(latParent,latChild, & 963 lonParent,lonChild, & 964 masksrc,matrix,src_add,dst_add) 965 966 END SELECT 967 ! 968 ENDIF 969 ! 970 SELECT CASE(TRIM(interp_type)) 971 CASE('bilinear') 972 CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, & 973 matrix,src_add,dst_add) 974 CASE('bicubic') 975 CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, & 976 matrix,src_add,dst_add) 977 END SELECT 978 ! 979 IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), & 980 G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 981 ! 982 IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1) 983 ! 984 dimnames(1)='x' 985 dimnames(2)='y' 986 dimnames(3)='time_counter' 987 ! 988 CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),& 989 Child_file,tabinterp3d,t,'float') 990 ! 991 DEALLOCATE(tabinterp3d) 992 DEALLOCATE(tabvar3d) 993 !end loop on time 994 END DO 995 ! 996 DEALLOCATE(detected_pts) 997 IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add) 998 DEALLOCATE( masksrc) 999 1000 ELSE 1001 1002 dimnames(1)='x' 1003 dimnames(2)='y' 1004 IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht' 1005 IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu' 1006 IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv' 1007 IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw' 1008 IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z' 1009 ! 1010 ! loop on vertical levels 1011 ! 1012 DO nb = 1,deptht 1013 ! 1014 ALLOCATE(masksrc(numlon,numlat)) 1015 ALLOCATE(detected_pts(numlon,numlat,N)) 1016 ! 1017 ! point detection et level n 1018 ! 1019 land_level = .FALSE. 1020 IF( Extrapolation ) THEN 1021 IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE. 1022 ENDIF 1023 1024 1025 IF ( land_level ) THEN 1026 ! 1027 WRITE(*,*) 'only land points on level ',nb 1028 ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 1029 tabinterp3d = 0.e0 1030 ! 1031 CALL Write_Ncdf_var(TRIM(varname),dimnames, & 1032 Child_file,tabinterp3d,nb,'float') 1033 DEALLOCATE(tabinterp3d) 1034 ! 1035 ELSE 1036 ! 1037 ALLOCATE(tabvar01(numlon,numlat,1)) ! level k 1038 IF(Extrapolation) ALLOCATE(tabvar02(numlon,numlat,1)) ! level k-1 1039 IF(Extrapolation) ALLOCATE(tabvar03(numlon,numlat,1)) ! level k-2 1040 ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 1041 ! 1042 IF(Extrapolation) THEN 1043 IF(nb==1) THEN 1044 CALL Read_Ncdf_var(varname,filename,tabvar01,nb) 1045 WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0. 1046 ELSE IF (nb==2) THEN 1047 CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1) 1048 CALL Read_Ncdf_var(varname,filename,tabvar01,nb) 1049 WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0. 1050 WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0. 1051 ELSE 1052 CALL Read_Ncdf_var(varname,filename,tabvar03,nb-2) 1053 CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1) 1054 CALL Read_Ncdf_var(varname,filename,tabvar01,nb) 1055 WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0. 1056 WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0. 1057 WHERE( tabvar03 .GE. 1.e+20 ) tabvar03 = 0. 1058 ENDIF 1059 ! 1060 CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb) 1061 1062 ALLOCATE(tabvar1(numlon,numlat,1,1)) ! level k 1063 ALLOCATE(tabvar2(numlon,numlat,1,1)) ! level k-1 1064 ALLOCATE(tabvar3(numlon,numlat,1,1)) ! level k-2 1065 tabvar1(:,:,1,1) = tabvar01(:,:,1) 1066 tabvar2(:,:,1,1) = tabvar02(:,:,1) 1067 tabvar3(:,:,1,1) = tabvar03(:,:,1) 1068 CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,& 1069 tabvar3,G0,depth,masksrc,nb,'posvar') 1070 tabvar01(:,:,1) = tabvar1(:,:,1,1) 1071 tabvar02(:,:,1) = tabvar2(:,:,1,1) 1072 tabvar03(:,:,1) = tabvar3(:,:,1,1) 1073 DEALLOCATE(tabvar1,tabvar2,tabvar3) 1074 DEALLOCATE(tabvar02,tabvar03) 1075 1076 ELSE 1077 CALL Read_Ncdf_var(varname,filename,tabvar01,nb) 1078 IF(MAXVAL(tabvar01(:,:,1))==0.) land_level = .TRUE. 1079 masksrc = .TRUE. 1080 ENDIF 1081 1082 IF( Extrapolation ) THEN 1083 WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for vertical level = ',nb 1084 ELSE 1085 WRITE(*,*) 'interpolation ',TRIM(varname),' for vertical level = ',nb 1086 ENDIF 1087 ! 1088 SELECT CASE(TRIM(interp_type)) 1089 CASE('bilinear') 1090 CALL get_remap_matrix(latParent,latChild, & 1091 lonParent,lonChild, & 1092 masksrc,matrix,src_add,dst_add) 1093 1094 CASE('bicubic') 1095 CALL get_remap_bicub(latParent,latChild, & 1096 lonParent,lonChild, & 1097 masksrc,matrix,src_add,dst_add) 1098 ! 1099 END SELECT 1100 ! 1101 ! 1102 SELECT CASE(TRIM(interp_type)) 1103 ! 1104 CASE('bilinear') 1105 CALL make_remap(tabvar01(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, & 1106 matrix,src_add,dst_add) 1107 CASE('bicubic') 1108 CALL make_bicubic_remap(tabvar01(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, & 1109 matrix,src_add,dst_add) 1110 END SELECT 1111 ! 1112 IF( conservation ) CALL Correctforconservation(tabvar01(:,:,1),tabinterp3d(:,:,1), & 1113 G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 1114 1115 ! 1116 ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 1117 tabinterp4d(:,:,1,1) = tabinterp3d(:,:,1) 1118 IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb) 1119 tabinterp3d(:,:,1) = tabinterp4d(:,:,1,1) 1120 DEALLOCATE(tabinterp4d) 1121 ! 1122 IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,nb) 1123 ! 1124 CALL Write_Ncdf_var(TRIM(varname),dimnames, & 1125 Child_file,tabinterp3d,nb,'float') 1126 ! 1127 DEALLOCATE(tabinterp3d) 1128 DEALLOCATE(tabvar01) 1129 ! 1130 ! 1131 ENDIF 1132 1133 ! 1134 IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add) 1135 DEALLOCATE( masksrc ) 1136 DEALLOCATE(detected_pts) 1137 ! 1138 ! end loop on vertical levels 1139 ! 1140 END DO 1141 1142 1143 ENDIF 887 1144 888 1145 CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file) -
utils/tools_r4.0-HEAD/NESTING/src/io_netcdf.f90
r10383 r14974 1604 1604 END FUNCTION Dims_Existence 1605 1605 ! 1606 LOGICAL FUNCTION DimUnlimited_Var( varname , filename ) 1607 ! 1608 CHARACTER(*),INTENT(in) :: varname,filename 1609 INTEGER :: ji, ncid, varid, status, numDims, unlimDimID 1610 INTEGER, DIMENSION(nf90_max_var_dims) :: VarDimIds 1611 ! 1612 status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid) 1613 IF (status/=nf90_noerr) THEN 1614 WRITE(*,*)"unable to open netcdf file : ",TRIM(filename) 1615 STOP 1616 ENDIF 1617 status = nf90_inquire(ncid, unlimiteddimid = unlimdimid) 1618 ! 1619 status = nf90_inq_varid(ncid,TRIM(varname),varid) 1620 ! 1621 status = nf90_inquire_variable(ncid, varid , ndims = numdims) 1622 status = nf90_inquire_variable(ncid, varid , dimids = vardimids(:numdims)) 1623 1624 DimUnlimited_Var = .FALSE. 1625 DO ji = 1, numdims 1626 IF( vardimids(ji) == unlimdimid ) DimUnlimited_Var = .TRUE. 1627 ENDDO 1628 1629 RETURN 1630 ! 1631 END FUNCTION DimUnlimited_Var 1632 1606 1633 END MODULE io_netcdf
Note: See TracChangeset
for help on using the changeset viewer.
