[Dart-dev] DART/branches Revision: 11130

[dart at ucar.edu]

mizzi at ucar.edu
2017-02-22 11:20:13 -0700 (Wed, 22 Feb 2017)
34
Update WRF-Chem/DART repository.




Modified: DART/branches/mizzi/observations/MOPITT_CO/mopitt_ascii_to_obs.f90
===================================================================
--- DART/branches/mizzi/observations/MOPITT_CO/mopitt_ascii_to_obs.f90	2017-02-22 18:19:01 UTC (rev 11129)
+++ DART/branches/mizzi/observations/MOPITT_CO/mopitt_ascii_to_obs.f90	2017-02-22 18:20:13 UTC (rev 11130)
@@ -174,6 +174,7 @@
 character*4             :: chr_year
 character*129           :: filedir, filename, copy_meta_data, filen
 character*129           :: transform_typ
+character*129           :: MOPITT_CO_retrieval_type
 !
 ! QOR/CPSR variables
 integer                                        :: info,nlvls_trc,qstatus
@@ -187,8 +188,15 @@
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
-namelist /create_mopitt_obs_nml/filedir,filename,year,month,day,hour,bin_beg, bin_end
+! MOPITT_CO_retrieval_type:
+!     RAWR - retrievals in VMR (ppb) units
+!     RETR - retrievals in log10(VMR ([ ])) units
+!     QOR  - quasi-optimal retrievals
+!     CPSR - compact phase space retrievals
 !
+namelist /create_mopitt_obs_nml/filedir,filename,year,month,day,hour,bin_beg, bin_end, &
+         MOPITT_CO_retrieval_type
+!
 ! Set constants
 ln_10=log(10.)
 pi=4.*atan(1.)
@@ -543,8 +551,8 @@
         xg_lat(i,j)=xg_lat(i,j)/xg_twt(i,j)
 !
 ! Skip for SINGLE_CLUSTER   
-        if((xg_lon(i,j).lt.-97. .or. xg_lon(i,j).gt.-93.) .or. &
-           (xg_lat(i,j).lt.38.  .or. xg_lat(i,j).gt.42.)) cycle
+!        if((xg_lon(i,j).lt.-97. .or. xg_lon(i,j).gt.-93.) .or. &
+!           (xg_lat(i,j).lt.38.  .or. xg_lat(i,j).gt.42.)) cycle
 !
         do k=1,mop_dim
            if(xg_norm(i,j,k).eq.0) cycle
@@ -642,151 +650,323 @@
 !             print *, 'cov ',k,(xg_ret_cov(i,j,k,l),l=kstr,mop_dim)
 !           enddo
 !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! QOR CODE BLOCK
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+        if(trim(MOPITT_CO_retrieval_type) .eq. 'QOR') then
+!
+! Calculate SVD of ret_cov (Z=U_xxx * SV_xxx * VT_xxx)
+           allocate(Z(nlvls,nlvls),SV_cov(nlvls),SV(nlvls,nlvls))
+           allocate(U_cov(nlvls,nlvls),UT_cov(nlvls,nlvls),V_cov(nlvls,nlvls),VT_cov(nlvls,nlvls))
+           allocate(rs_avg_k(nlvls,nlvls),rs_cov(nlvls,nlvls),rs_x_r(nlvls),rs_x_p(nlvls))       
+           allocate(rr_avg_k(nlvls,nlvls),rr_cov(nlvls,nlvls),rr_x_r(nlvls),rr_x_p(nlvls))       
+           allocate(ZL(nlvls,nlvls),ZR(nlvls,nlvls),ZV(nlvls))
+           Z(1:nlvls,1:nlvls)=dble(xg_ret_cov(i,j,kstr:mop_dim,kstr:mop_dim))
+           call dgesvd('A','A',nlvls,nlvls,Z,nlvls,SV_cov,U_cov,nlvls,VT_cov,nlvls,wrk,lwrk,info)
+           nlvls_trc=0
+           do k=1,nlvls
+              if(SV_cov(k).ge.eps_tol) then
+                 nlvls_trc=k
+              else
+                 SV_cov(k)=0
+                 U_cov(:,k)=0. 
+                 VT_cov(k,:)=0.
+              endif 
+           enddo
+!              print *,'nlvls_trc ',nlvls_trc
+!              print *, 'SV ',SV_cov(:)
+!
+! Scale the singular vectors (NO SCALE/SCALE)     
+           do k=1,nlvls_trc
+              U_cov(:,k)=U_cov(:,k)/sqrt(SV_cov(k))
+           enddo
+!              print *, 'nlvls_trc ',nlvls_trc
+!              print *, 'SV ',SV_cov(:)
+!          
+           call mat_transpose(U_cov,UT_cov,nlvls,nlvls)
+           call mat_transpose(VT_cov,V_cov,nlvls,nlvls)
+           call vec_to_mat(SV_cov,SV,nlvls)
+!              do k=1,nlvls
+!                print *, 'U ',k,(U_cov(k,l),l=1,nlvls)
+!              enddo
+!              do k=1,nlvls
+!                print *, 'UT ',k,(UT_cov(k,l),l=1,nlvls)
+!              enddo
+!
+! Rotate terms in the forward operator
+           ZL(1:nlvls,1:nlvls)=dble(xg_avg_k(i,j,kstr:mop_dim,kstr:mop_dim))
+           call mat_prd(UT_cov(1:nlvls,1:nlvls),ZL(1:nlvls,1:nlvls), &
+           rs_avg_k(1:nlvls,1:nlvls),nlvls,nlvls,nlvls,nlvls)
+!              do k=1,nlvls
+!                print *, 'UT ',k,(UT_cov(k,l),l=1,nlvls)
+!              enddo
+!              do k=kstr,mop_dim
+!                print *, 'xg_avg_k ',k,(xg_avg_k(i,j,k,l),l=kstr,mop_dim)
+!              enddo
+!              do k=1,nlvls