Changeset 9211 for branches/UKMO
- Timestamp:
- 2018-01-11T16:43:41+01:00 (6 years ago)
- Location:
- branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS
- Files:
-
- 2 edited
-
diaobs.F90 (modified) (9 diffs)
-
obs_oper.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r9205 r9211 747 747 & kdailyavtypes = nn_profdavtypes ) 748 748 749 ! Is allocating and deallocating repeatedly in a loop good practice?750 749 DEALLOCATE( llvar ) 751 750 CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zglam ) … … 889 888 INTEGER :: ji, jj ! Loop counters 890 889 REAL(wp) :: tiny ! small number 891 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 892 & zprofvar1, & ! Model values for 1st variable in a prof ob 893 & zprofvar2 ! Model values for 2nd variable in a prof ob 894 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 895 & zprofmask1, & ! Mask associated with zprofvar1 896 & zprofmask2 ! Mask associated with zprofvar2 890 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 891 & zprofvar ! Model values for variables in a prof ob 892 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 893 & zprofmask ! Mask associated with zprofvar 897 894 REAL(wp), POINTER, DIMENSION(:,:) :: & 898 895 & zsurfvar, & ! Model values equivalent to surface ob. 899 896 & zsurfmask ! Mask associated with surface variable 900 REAL(wp), POINTER, DIMENSION(:,:) :: & 901 & zglam1, & ! Model longitudes for prof variable 1 902 & zglam2, & ! Model longitudes for prof variable 2 903 & zgphi1, & ! Model latitudes for prof variable 1 904 & zgphi2 ! Model latitudes for prof variable 2 897 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 898 & zglam, & ! Model longitudes for prof variables 899 & zgphi ! Model latitudes for prof variables 905 900 LOGICAL :: llog10 ! Perform log10 transform of variable 906 901 907 908 !Allocate local work arrays909 CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )910 CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )911 CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )912 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )913 CALL wrk_alloc( jpi, jpj, zsurfvar )914 CALL wrk_alloc( jpi, jpj, zsurfmask )915 CALL wrk_alloc( jpi, jpj, zglam1 )916 CALL wrk_alloc( jpi, jpj, zglam2 )917 CALL wrk_alloc( jpi, jpj, zgphi1 )918 CALL wrk_alloc( jpi, jpj, zgphi2 )919 902 920 903 IF(lwp) THEN … … 935 918 DO jtype = 1, nproftypes 936 919 920 ! Allocate local work arrays 921 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 922 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 923 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 924 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 925 926 ! Defaults which might change 927 DO jvar = 1, profdataqc(jtype)%nvar 928 zprofmask(:,:,:,jvar) = tmask(:,:,:) 929 zglam(:,:,jvar) = glamt(:,:) 930 zgphi(:,:,jvar) = gphit(:,:) 931 END DO 932 937 933 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 934 938 935 CASE('prof') 939 zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 940 zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 941 zprofmask1(:,:,:) = tmask(:,:,:) 942 zprofmask2(:,:,:) = tmask(:,:,:) 943 zglam1(:,:) = glamt(:,:) 944 zglam2(:,:) = glamt(:,:) 945 zgphi1(:,:) = gphit(:,:) 946 zgphi2(:,:) = gphit(:,:) 936 zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 937 zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 938 947 939 CASE('vel') 948 zprofvar 1(:,:,:) = un(:,:,:)949 zprofvar 2(:,:,:) = vn(:,:,:)950 zprofmask 1(:,:,:) = umask(:,:,:)951 zprofmask 2(:,:,:) = vmask(:,:,:)952 zglam 1(:,:) = glamu(:,:)953 zglam 2(:,:) = glamv(:,:)954 zgphi 1(:,:) = gphiu(:,:)955 zgphi 2(:,:) = gphiv(:,:)940 zprofvar(:,:,:,1) = un(:,:,:) 941 zprofvar(:,:,:,2) = vn(:,:,:) 942 zprofmask(:,:,:,1) = umask(:,:,:) 943 zprofmask(:,:,:,2) = vmask(:,:,:) 944 zglam(:,:,1) = glamu(:,:) 945 zglam(:,:,2) = glamv(:,:) 946 zgphi(:,:,1) = gphiu(:,:) 947 zgphi(:,:,2) = gphiv(:,:) 956 948 957 949 CASE('plchltot') 958 950 #if defined key_hadocc 959 951 ! Chlorophyll from HadOCC 960 zprofvar 1(:,:,:) = HADOCC_CHL(:,:,:)952 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 961 953 #elif defined key_medusa && defined key_foam_medusa 962 954 ! Add non-diatom and diatom chlorophyll from MEDUSA 963 zprofvar 1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)955 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 964 956 #elif defined key_fabm 965 957 ! Add all chlorophyll groups from ERSEM 966 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &967 & trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)958 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + & 959 & trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4) 968 960 #else 969 961 CALL ctl_stop( ' Trying to run plchltot observation operator', & 970 962 & ' but no biogeochemical model appears to have been defined' ) 971 963 #endif 972 zprofmask1(:,:,:) = tmask(:,:,:)973 964 ! Take the log10 where we can, otherwise exclude 974 965 tiny = 1.0e-20 975 WHERE(zprofvar 1(:,:,:) > tiny .AND. zprofvar1(:,:,:) /= obfillflt )976 zprofvar 1(:,:,:) = LOG10(zprofvar1(:,:,:))966 WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 967 zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:)) 977 968 ELSEWHERE 978 zprofvar 1(:,:,:) = obfillflt979 zprofmask 1(:,:,:) = 0969 zprofvar(:,:,:,:) = obfillflt 970 zprofmask(:,:,:,:) = 0 980 971 END WHERE 981 zglam1(:,:) = glamt(:,:)982 zgphi1(:,:) = gphit(:,:)983 972 984 973 CASE('pchltot') 985 974 #if defined key_hadocc 986 975 ! Chlorophyll from HadOCC 987 zprofvar 1(:,:,:) = HADOCC_CHL(:,:,:)976 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 988 977 #elif defined key_medusa && defined key_foam_medusa 989 978 ! Add non-diatom and diatom chlorophyll from MEDUSA 990 zprofvar 1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)979 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 991 980 #elif defined key_fabm 992 981 ! Add all chlorophyll groups from ERSEM 993 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &994 & trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)982 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + & 983 & trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4) 995 984 #else 996 985 CALL ctl_stop( ' Trying to run pchltot observation operator', & 997 986 & ' but no biogeochemical model appears to have been defined' ) 998 987 #endif 999 zprofmask1(:,:,:) = tmask(:,:,:)1000 zglam1(:,:) = glamt(:,:)1001 zgphi1(:,:) = gphit(:,:)1002 988 1003 989 CASE('pno3') 1004 990 #if defined key_hadocc 1005 991 ! Dissolved inorganic nitrogen from HadOCC 1006 zprofvar 1(:,:,:) = trn(:,:,:,jp_had_nut)992 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 1007 993 #elif defined key_medusa && defined key_foam_medusa 1008 994 ! Dissolved inorganic nitrogen from MEDUSA 1009 zprofvar 1(:,:,:) = trn(:,:,:,jpdin)995 zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 1010 996 #elif defined key_fabm 1011 997 ! Nitrate from ERSEM 1012 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_n3n)998 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n3n) 1013 999 #else 1014 1000 CALL ctl_stop( ' Trying to run pno3 observation operator', & 1015 1001 & ' but no biogeochemical model appears to have been defined' ) 1016 1002 #endif 1017 zprofmask1(:,:,:) = tmask(:,:,:)1018 zglam1(:,:) = glamt(:,:)1019 zgphi1(:,:) = gphit(:,:)1020 1003 1021 1004 CASE('psi4') … … 1025 1008 #elif defined key_medusa && defined key_foam_medusa 1026 1009 ! Silicate from MEDUSA 1027 zprofvar 1(:,:,:) = trn(:,:,:,jpsil)1010 zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 1028 1011 #elif defined key_fabm 1029 1012 ! Silicate from ERSEM 1030 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_n5s)1013 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n5s) 1031 1014 #else 1032 1015 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1033 1016 & ' but no biogeochemical model appears to have been defined' ) 1034 1017 #endif 1035 zprofmask1(:,:,:) = tmask(:,:,:)1036 zglam1(:,:) = glamt(:,:)1037 zgphi1(:,:) = gphit(:,:)1038 1018 1039 1019 CASE('ppo4') … … 1046 1026 #elif defined key_fabm 1047 1027 ! Phosphate from ERSEM 1048 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_n1p)1028 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n1p) 1049 1029 #else 1050 1030 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1051 1031 & ' but no biogeochemical model appears to have been defined' ) 1052 1032 #endif 1053 zprofmask1(:,:,:) = tmask(:,:,:)1054 zglam1(:,:) = glamt(:,:)1055 zgphi1(:,:) = gphit(:,:)1056 1033 1057 1034 CASE('pdic') 1058 1035 #if defined key_hadocc 1059 1036 ! Dissolved inorganic carbon from HadOCC 1060 zprofvar 1(:,:,:) = trn(:,:,:,jp_had_dic)1037 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 1061 1038 #elif defined key_medusa && defined key_foam_medusa 1062 1039 ! Dissolved inorganic carbon from MEDUSA 1063 zprofvar 1(:,:,:) = trn(:,:,:,jpdic)1040 zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 1064 1041 #elif defined key_fabm 1065 1042 ! Dissolved inorganic carbon from ERSEM 1066 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_o3c)1043 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3c) 1067 1044 #else 1068 1045 CALL ctl_stop( ' Trying to run pdic observation operator', & 1069 1046 & ' but no biogeochemical model appears to have been defined' ) 1070 1047 #endif 1071 zprofmask1(:,:,:) = tmask(:,:,:)1072 zglam1(:,:) = glamt(:,:)1073 zgphi1(:,:) = gphit(:,:)1074 1048 1075 1049 CASE('palk') 1076 1050 #if defined key_hadocc 1077 1051 ! Alkalinity from HadOCC 1078 zprofvar 1(:,:,:) = trn(:,:,:,jp_had_alk)1052 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 1079 1053 #elif defined key_medusa && defined key_foam_medusa 1080 1054 ! Alkalinity from MEDUSA 1081 zprofvar 1(:,:,:) = trn(:,:,:,jpalk)1055 zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 1082 1056 #elif defined key_fabm 1083 1057 ! Alkalinity from ERSEM 1084 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_o3a)1058 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3a) 1085 1059 #else 1086 1060 CALL ctl_stop( ' Trying to run palk observation operator', & 1087 1061 & ' but no biogeochemical model appears to have been defined' ) 1088 1062 #endif 1089 zprofmask1(:,:,:) = tmask(:,:,:)1090 zglam1(:,:) = glamt(:,:)1091 zgphi1(:,:) = gphit(:,:)1092 1063 1093 1064 CASE('pph') … … 1097 1068 #elif defined key_medusa && defined key_foam_medusa 1098 1069 ! pH from MEDUSA 1099 zprofvar 1(:,:,:) = f3_pH(:,:,:)1070 zprofvar(:,:,:,1) = f3_pH(:,:,:) 1100 1071 #elif defined key_fabm 1101 1072 ! pH from ERSEM 1102 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_o3ph)1073 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3ph) 1103 1074 #else 1104 1075 CALL ctl_stop( ' Trying to run pph observation operator', & 1105 1076 & ' but no biogeochemical model appears to have been defined' ) 1106 1077 #endif 1107 zprofmask1(:,:,:) = tmask(:,:,:)1108 zglam1(:,:) = glamt(:,:)1109 zgphi1(:,:) = gphit(:,:)1110 1078 1111 1079 CASE('po2') … … 1115 1083 #elif defined key_medusa && defined key_foam_medusa 1116 1084 ! Oxygen from MEDUSA 1117 zprofvar 1(:,:,:) = trn(:,:,:,jpoxy)1085 zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 1118 1086 #elif defined key_fabm 1119 1087 ! Oxygen from ERSEM 1120 zprofvar 1(:,:,:) = trn(:,:,:,jp_fabm_o2o)1088 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o2o) 1121 1089 #else 1122 1090 CALL ctl_stop( ' Trying to run po2 observation operator', & 1123 1091 & ' but no biogeochemical model appears to have been defined' ) 1124 1092 #endif 1125 zprofmask1(:,:,:) = tmask(:,:,:)1126 zglam1(:,:) = glamt(:,:)1127 zgphi1(:,:) = gphit(:,:)1128 1093 1129 1094 CASE DEFAULT 1130 1095 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 1096 1131 1097 END SELECT 1132 1133 IF ( ( TRIM(cobstypesprof(jtype)) /= 'prof' ) .AND. ( TRIM(cobstypesprof(jtype)) /= 'vel' ) ) THEN1134 zprofvar2(:,:,:) = zprofvar1(:,:,:)1135 zprofmask2(:,:,:) = zprofmask1(:,:,:)1136 zglam2(:,:) = zglam1(:,:)1137 zgphi2(:,:) = zgphi1(:,:)1138 ENDIF1139 1140 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,&1141 & nit000, idaystp, &1142 & zprofvar1, zprofvar2, &1143 & fsdept(:,:,:), fsdepw(:,:,:), & 1144 & zprofmask1, zprofmask2, &1145 & zglam1, zglam2, zgphi1, zgphi2, &1146 & nn_1dint, nn_2dint_default, &1147 & kdailyavtypes = nn_profdavtypes)1098 1099 DO jvar = 1, profdataqc(jtype)%nvar 1100 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 1101 & nit000, idaystp, jvar, & 1102 & zprofvar(:,:,:,jvar), & 1103 & fsdept(:,:,:), fsdepw(:,:,:), & 1104 & zprofmask(:,:,:,jvar), & 1105 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1106 & nn_1dint, nn_2dint_default, & 1107 & kdailyavtypes = nn_profdavtypes ) 1108 END DO 1109 1110 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1111 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1112 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1113 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 1148 1114 1149 1115 END DO … … 1152 1118 1153 1119 IF ( nsurftypes > 0 ) THEN 1120 1121 !Allocate local work arrays 1122 CALL wrk_alloc( jpi, jpj, zsurfvar ) 1123 CALL wrk_alloc( jpi, jpj, zsurfmask ) 1154 1124 1155 1125 DO jtype = 1, nsurftypes … … 1414 1384 END DO 1415 1385 1386 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 1387 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 1388 1416 1389 ENDIF 1417 1418 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )1419 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )1420 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )1421 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )1422 CALL wrk_dealloc( jpi, jpj, zsurfvar )1423 CALL wrk_dealloc( jpi, jpj, zsurfmask )1424 CALL wrk_dealloc( jpi, jpj, zglam1 )1425 CALL wrk_dealloc( jpi, jpj, zglam2 )1426 CALL wrk_dealloc( jpi, jpj, zgphi1 )1427 CALL wrk_dealloc( jpi, jpj, zgphi2 )1428 1390 1429 1391 END SUBROUTINE dia_obs -
branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r9186 r9211 36 36 & sbc_dcy, nday_qsr 37 37 USE obs_grid, ONLY : & 38 & obs_level_search 38 & obs_level_search 39 USE dom_oce ! Ocean space and time domain variables 39 40 40 41 IMPLICIT NONE … … 60 61 61 62 62 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, 63 & kit000, kdaystp, 64 & pvar 1, pvar2, pgdept, pgdepw,&65 & pmask 1, pmask2, &66 & plam 1, plam2, pphi1, pphi2,&63 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 64 & kit000, kdaystp, kvar, & 65 & pvar, pgdept, pgdepw, & 66 & pmask, & 67 & plam, pphi, & 67 68 & k1dint, k2dint, kdailyavtypes ) 68 69 … … 134 135 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 135 136 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 137 INTEGER, INTENT(IN) :: kvar ! Number of variable in prodatqc 136 138 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 137 & pvar1, & ! Model field 1 138 & pvar2, & ! Model field 2 139 & pmask1, & ! Land-sea mask 1 140 & pmask2 ! Land-sea mask 2 139 & pvar, & ! Model field for variable 140 & pmask ! Land-sea mask for variable 141 141 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 142 & plam1, & ! Model longitudes for variable 1 143 & plam2, & ! Model longitudes for variable 2 144 & pphi1, & ! Model latitudes for variable 1 145 & pphi2 ! Model latitudes for variable 2 142 & plam, & ! Model longitudes for variable 143 & pphi ! Model latitudes for variable 146 144 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 147 145 & pgdept, & ! Model array of depth T levels … … 166 164 & idailyavtypes 167 165 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 168 & igrdi1, & 169 & igrdi2, & 170 & igrdj1, & 171 & igrdj2 166 & igrdi, & 167 & igrdj 172 168 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 173 169 … … 176 172 REAL(KIND=wp) :: zdaystp 177 173 REAL(KIND=wp), DIMENSION(kpk) :: & 178 & zobsmask1, &179 & zobsmask2, &180 174 & zobsk, & 181 175 & zobs2k 182 176 REAL(KIND=wp), DIMENSION(2,2,1) :: & 183 177 & zweig1, & 184 & zweig2, &185 178 & zweig 186 179 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 187 & zmask1, & 188 & zmask2, & 189 & zint1, & 190 & zint2, & 191 & zinm1, & 192 & zinm2, & 180 & zmask, & 181 & zint, & 182 & zinm, & 193 183 & zgdept, & 194 184 & zgdepw 195 185 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 196 & zglam1, & 197 & zglam2, & 198 & zgphi1, & 199 & zgphi2 200 REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2 186 & zglam, & 187 & zgphi 188 REAL(KIND=wp), DIMENSION(1) :: zmsk 201 189 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 202 190 … … 230 218 DO jj = 1, jpj 231 219 DO ji = 1, jpi 232 prodatqc%vdmean(ji,jj,jk,1) = 0.0 233 IF ( prodatqc%nvar == 2 ) THEN 234 prodatqc%vdmean(ji,jj,jk,2) = 0.0 235 ENDIF 220 prodatqc%vdmean(ji,jj,jk,kvar) = 0.0 236 221 END DO 237 222 END DO … … 242 227 DO jj = 1, jpj 243 228 DO ji = 1, jpi 244 ! Increment field 1 for computing daily mean 245 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 246 & + pvar1(ji,jj,jk) 247 IF ( prodatqc%nvar == 2 ) THEN 248 ! Increment field 2 for computing daily mean 249 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 250 & + pvar2(ji,jj,jk) 251 ENDIF 229 ! Increment field for computing daily mean 230 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 231 & + pvar(ji,jj,jk) 252 232 END DO 253 233 END DO … … 262 242 DO jj = 1, jpj 263 243 DO ji = 1, jpi 264 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 265 & * zdaystp 266 IF ( prodatqc%nvar == 2 ) THEN 267 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 268 & * zdaystp 269 ENDIF 244 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 245 & * zdaystp 270 246 END DO 271 247 END DO … … 277 253 ! Get the data for interpolation 278 254 ALLOCATE( & 279 & igrdi1(2,2,ipro), & 280 & igrdi2(2,2,ipro), & 281 & igrdj1(2,2,ipro), & 282 & igrdj2(2,2,ipro), & 283 & zglam1(2,2,ipro), & 284 & zglam2(2,2,ipro), & 285 & zgphi1(2,2,ipro), & 286 & zgphi2(2,2,ipro), & 287 & zmask1(2,2,kpk,ipro), & 288 & zmask2(2,2,kpk,ipro), & 289 & zint1(2,2,kpk,ipro), & 290 & zint2(2,2,kpk,ipro), & 255 & igrdi(2,2,ipro), & 256 & igrdj(2,2,ipro), & 257 & zglam(2,2,ipro), & 258 & zgphi(2,2,ipro), & 259 & zmask(2,2,kpk,ipro), & 260 & zint(2,2,kpk,ipro), & 291 261 & zgdept(2,2,kpk,ipro), & 292 262 & zgdepw(2,2,kpk,ipro) & … … 295 265 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 296 266 iobs = jobs - prodatqc%nprofup 297 igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 298 igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 299 igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 300 igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 301 igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 302 igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 303 igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 304 igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 305 IF ( prodatqc%nvar == 2 ) THEN 306 igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 307 igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 308 igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 309 igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 310 igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 311 igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 312 igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 313 igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 314 ENDIF 267 igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1 268 igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1 269 igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1 270 igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar) 271 igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar) 272 igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1 273 igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar) 274 igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar) 315 275 END DO 316 276 … … 319 279 zgdepw(:,:,:,:) = 0.0 320 280 321 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 322 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 323 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 324 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 ) 325 326 IF ( prodatqc%nvar == 2 ) THEN 327 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 328 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 329 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 330 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 331 ENDIF 332 333 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 334 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 281 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam ) 282 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 283 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar, zint ) 284 285 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 286 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 335 287 336 288 ! At the end of the day also get interpolated means 337 289 IF ( ld_dailyav .AND. idayend == 0 ) THEN 338 290 339 ALLOCATE( & 340 & zinm1(2,2,kpk,ipro), & 341 & zinm2(2,2,kpk,ipro) & 342 & ) 343 344 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 345 & prodatqc%vdmean(:,:,:,1), zinm1 ) 346 IF ( prodatqc%nvar == 2 ) THEN 347 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 348 & prodatqc%vdmean(:,:,:,2), zinm2 ) 349 ENDIF 291 ALLOCATE( zinm(2,2,kpk,ipro) ) 292 293 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 294 & prodatqc%vdmean(:,:,:,kvar), zinm ) 350 295 351 296 ENDIF … … 382 327 ! Horizontal weights 383 328 ! Masked values are calculated later. 384 IF ( prodatqc%npvend(jobs, 1) > 0 ) THEN329 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 385 330 386 331 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 387 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 388 & zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 389 390 ENDIF 391 392 IF ( prodatqc%nvar == 2 ) THEN 393 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 394 395 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 396 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 397 & zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 398 399 ENDIF 400 ENDIF 401 402 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 332 & zglam(:,:,iobs), zgphi(:,:,iobs), & 333 & zmask(:,:,1,iobs), zweig1, zmsk ) 334 335 ENDIF 336 337 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 403 338 404 339 zobsk(:) = obfillflt … … 420 355 IF ( k1dint == 1 ) THEN 421 356 CALL obs_int_z1d_spl( kpk, & 422 & zinm 1(iin,ijn,:,iobs), &357 & zinm(iin,ijn,:,iobs), & 423 358 & zobs2k, zgdept(iin,ijn,:,iobs), & 424 & zmask 1(iin,ijn,:,iobs))359 & zmask(iin,ijn,:,iobs)) 425 360 ENDIF 426 361 427 362 CALL obs_level_search(kpk, & 428 363 & zgdept(iin,ijn,:,iobs), & 429 & inum_obs, prodatqc%var( 1)%vdep(ista:iend), &364 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 430 365 & iv_indic) 431 366 432 367 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 433 & prodatqc%var( 1)%vdep(ista:iend), &434 & zinm 1(iin,ijn,:,iobs), &368 & prodatqc%var(kvar)%vdep(ista:iend), & 369 & zinm(iin,ijn,:,iobs), & 435 370 & zobs2k, interp_corner(iin,ijn,:), & 436 371 & zgdept(iin,ijn,:,iobs), & 437 & zmask 1(iin,ijn,:,iobs))372 & zmask(iin,ijn,:,iobs)) 438 373 439 374 ENDDO … … 447 382 448 383 ! vertically interpolate all 4 corners 449 ista = prodatqc%npvsta(jobs, 1)450 iend = prodatqc%npvend(jobs, 1)384 ista = prodatqc%npvsta(jobs,kvar) 385 iend = prodatqc%npvend(jobs,kvar) 451 386 inum_obs = iend - ista + 1 452 387 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) … … 456 391 IF ( k1dint == 1 ) THEN 457 392 CALL obs_int_z1d_spl( kpk, & 458 & zint 1(iin,ijn,:,iobs),&393 & zint(iin,ijn,:,iobs),& 459 394 & zobs2k, zgdept(iin,ijn,:,iobs), & 460 & zmask 1(iin,ijn,:,iobs))395 & zmask(iin,ijn,:,iobs)) 461 396 462 397 ENDIF … … 464 399 CALL obs_level_search(kpk, & 465 400 & zgdept(iin,ijn,:,iobs),& 466 & inum_obs, prodatqc%var( 1)%vdep(ista:iend), &401 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 467 402 & iv_indic) 468 403 469 404 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 470 & prodatqc%var( 1)%vdep(ista:iend), &471 & zint 1(iin,ijn,:,iobs), &405 & prodatqc%var(kvar)%vdep(ista:iend), & 406 & zint(iin,ijn,:,iobs), & 472 407 & zobs2k,interp_corner(iin,ijn,:), & 473 408 & zgdept(iin,ijn,:,iobs), & 474 & zmask 1(iin,ijn,:,iobs) )409 & zmask(iin,ijn,:,iobs) ) 475 410 476 411 ENDDO … … 496 431 DO ijn=1,2 497 432 498 depth_loop 1: DO ik=kpk,2,-1499 IF(zmask 1(iin,ijn,ik-1,iobs ) > 0.9 )THEN433 depth_loop: DO ik=kpk,2,-1 434 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 500 435 501 436 zweig(iin,ijn,1) = & 502 437 & zweig1(iin,ijn,1) * & 503 438 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 504 & - prodatqc%var( 1)%vdep(iend)),0._wp)439 & - prodatqc%var(kvar)%vdep(iend)),0._wp) 505 440 506 EXIT depth_loop 1441 EXIT depth_loop 507 442 508 443 ENDIF 509 444 510 ENDDO depth_loop 1445 ENDDO depth_loop 511 446 512 447 ENDDO … … 514 449 515 450 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 516 & prodatqc%var( 1)%vmod(iend:iend) )451 & prodatqc%var(kvar)%vmod(iend:iend) ) 517 452 518 453 ! Set QC flag for any observations found below the bottom 519 454 ! needed as the check here is more strict than that in obs_prep 520 IF (sum(zweig) == 0.0_wp) prodatqc%var( 1)%nvqc(iend:iend)=4455 IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4 521 456 522 457 ENDDO … … 524 459 DEALLOCATE(interp_corner,iv_indic) 525 460 526 ENDIF 527 528 IF ( prodatqc%nvar == 2 ) THEN 529 ! For the second variable 530 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 531 532 zobsk(:) = obfillflt 533 534 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 535 536 IF ( idayend == 0 ) THEN 537 ! Daily averaged data 538 539 ! vertically interpolate all 4 corners 540 ista = prodatqc%npvsta(jobs,2) 541 iend = prodatqc%npvend(jobs,2) 542 inum_obs = iend - ista + 1 543 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 544 545 DO iin=1,2 546 DO ijn=1,2 547 548 IF ( k1dint == 1 ) THEN 549 CALL obs_int_z1d_spl( kpk, & 550 & zinm2(iin,ijn,:,iobs), & 551 & zobs2k, zgdept(iin,ijn,:,iobs), & 552 & zmask2(iin,ijn,:,iobs)) 553 ENDIF 554 555 CALL obs_level_search(kpk, & 556 & zgdept(iin,ijn,:,iobs), & 557 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 558 & iv_indic) 559 560 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 561 & prodatqc%var(2)%vdep(ista:iend), & 562 & zinm2(iin,ijn,:,iobs), & 563 & zobs2k, interp_corner(iin,ijn,:), & 564 & zgdept(iin,ijn,:,iobs), & 565 & zmask2(iin,ijn,:,iobs)) 566 567 ENDDO 568 ENDDO 569 570 ENDIF !idayend 571 572 ELSE 573 574 ! Point data 575 576 ! vertically interpolate all 4 corners 577 ista = prodatqc%npvsta(jobs,2) 578 iend = prodatqc%npvend(jobs,2) 579 inum_obs = iend - ista + 1 580 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 581 DO iin=1,2 582 DO ijn=1,2 583 584 IF ( k1dint == 1 ) THEN 585 CALL obs_int_z1d_spl( kpk, & 586 & zint2(iin,ijn,:,iobs),& 587 & zobs2k, zgdept(iin,ijn,:,iobs), & 588 & zmask2(iin,ijn,:,iobs)) 589 590 ENDIF 591 592 CALL obs_level_search(kpk, & 593 & zgdept(iin,ijn,:,iobs),& 594 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 595 & iv_indic) 596 597 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 598 & prodatqc%var(2)%vdep(ista:iend), & 599 & zint2(iin,ijn,:,iobs), & 600 & zobs2k,interp_corner(iin,ijn,:), & 601 & zgdept(iin,ijn,:,iobs), & 602 & zmask2(iin,ijn,:,iobs) ) 603 604 ENDDO 605 ENDDO 606 607 ENDIF 608 609 !------------------------------------------------------------- 610 ! Compute the horizontal interpolation for every profile level 611 !------------------------------------------------------------- 612 613 DO ikn=1,inum_obs 614 iend=ista+ikn-1 615 616 zweig(:,:,1) = 0._wp 617 618 ! This code forces the horizontal weights to be 619 ! zero IF the observation is below the bottom of the 620 ! corners of the interpolation nodes, Or if it is in 621 ! the mask. This is important for observations near 622 ! steep bathymetry 623 DO iin=1,2 624 DO ijn=1,2 625 626 depth_loop2: DO ik=kpk,2,-1 627 IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN 628 629 zweig(iin,ijn,1) = & 630 & zweig2(iin,ijn,1) * & 631 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 632 & - prodatqc%var(2)%vdep(iend)),0._wp) 633 634 EXIT depth_loop2 635 636 ENDIF 637 638 ENDDO depth_loop2 639 640 ENDDO 641 ENDDO 642 643 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 644 & prodatqc%var(2)%vmod(iend:iend) ) 645 646 ! Set QC flag for any observations found below the bottom 647 ! needed as the check here is more strict than that in obs_prep 648 IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 649 650 ENDDO 651 652 DEALLOCATE(interp_corner,iv_indic) 653 654 ENDIF 655 ENDIF 461 ENDIF 656 462 657 463 ENDDO 658 464 659 465 ! Deallocate the data for interpolation 660 DEALLOCATE( & 661 & igrdi1, & 662 & igrdi2, & 663 & igrdj1, & 664 & igrdj2, & 665 & zglam1, & 666 & zglam2, & 667 & zgphi1, & 668 & zgphi2, & 669 & zmask1, & 670 & zmask2, & 671 & zint1, & 672 & zint2, & 466 DEALLOCATE( & 467 & igrdi, & 468 & igrdj, & 469 & zglam, & 470 & zgphi, & 471 & zmask, & 472 & zint, & 673 473 & zgdept, & 674 474 & zgdepw & … … 677 477 ! At the end of the day also get interpolated means 678 478 IF ( ld_dailyav .AND. idayend == 0 ) THEN 679 DEALLOCATE( & 680 & zinm1, & 681 & zinm2 & 682 & ) 479 DEALLOCATE( zinm ) 683 480 ENDIF 684 481 685 prodatqc%nprofup = prodatqc%nprofup + ipro 482 IF ( kvar == prodatqc%nvar ) THEN 483 prodatqc%nprofup = prodatqc%nprofup + ipro 484 ENDIF 686 485 687 486 END SUBROUTINE obs_prof_opt
Note: See TracChangeset
for help on using the changeset viewer.
