Last push before Zenobe dies

EnKF-MPI-TOPAZ/m_prep_4_EnKF.F90

@@ -288,7 +288,8 @@ contains
              write(cmem,'(i3.3)') iens
              tlevel = 1
              call get_mod_fld_new(trim('forecast'//cmem), readfld, iens,&
-                  'temp', 1, tlevel, nx, ny)
+                  'tn', 1, tlevel, nx, ny)
+             PRINT *, "FRANCOIS"
              if (tlevel == -1) then
                 if (master) then
                    print *, 'ERROR: get_mod_fld_new(): failed for "SST"'

conversion_uf/p_prep_obs_ORCA1.F90

@@ -39,6 +39,12 @@ program prep_obs
   ! Data vars
   integer, parameter           :: nx=362, ny=292 ! Unfortunately, still hard coded.
   real, dimension(nx,ny)       :: lons, lats, obsfld, errorfld, obsfldu, obsfldv, errorfldu, errorfldv
+
+  REAL                         :: obsmin, obsmax, errmin, errmax
+  REAL                         :: latmin_NH = 40.0
+  REAL                         :: latmax_NH = 90.0
+  REAL                         :: latmin_SH = 40.0 ! Same for SH (sign will be added)
+  REAL                         :: latmax_SH = 90.0
   integer, dimension(nx,ny)    :: mask
   integer                      :: obscnt       ! Counts nr of obs's.
   ! Other vars
@@ -117,15 +123,38 @@ program prep_obs
   read(dummy,'(A)') varname
   call getarg(3, dummy)
   read(dummy,*) rdate
-  if (      trim(varname) .ne.'rfb' .and. (trim(varname) .ne. 'vt_i') ) then
-    write(*,*) ' (prep_obs) Currently only the radar freeboard'
-    write(*,*) '            rfb    can be processed                      '
-    PRINT *, varname
-    stop " (prep_obs) Aborted"
-  endif
+
+
+  IF     ( TRIM(varname) .eq. 'rfb' ) THEN
+    WRITE(*,*) "Handling ", TRIM(varname)
+    ! Min and max values for error used to screen the data (any data with
+    ! standard deviation in between those values will be selected
+    obsmin = 0.01
+    obsmax = 10.0
+    errmin = 0.01
+    errmax = 1.0
+  ELSEIF ( TRIM(varname) .eq. 'vt_i') THEN
+    WRITE(*,*) "Handling ", TRIM(varname)
+    obsmin = 0.01
+    obsmax = 50.0
+    errmin = 0.01
+    errmax = 1.0
+  ELSEIF ( TRIM(varname) .eq. 'at_i') THEN
+    WRITE(*,*) "Handling ", TRIM(varname)
+    obsmin = 0.001
+    obsmax = 1.0
+    errmin = 0.001
+    errmax = 0.5
+  ELSE
+    WRITE(*,*) " (prep_obs) Currently only the variables rfb  (sea ice freeboard),"
+    WRITE(*,*) "                                         at_i (total sea ice volume)"
+    WRITE(*,*) "                                         at_i (total sea ice concentration)"
+    WRITE(*,*) "            can be processed                      "
+    STOP       "(prep_obs) Aborted"
+  ENDIF
 
   ! User info
-  write(*,*) ' (prep_obs) Extracting '//trim(varname)//' from '//trim(filename)
+  WRITE(*,*)  "(prep_obs) Extracting "//TRIM(varname)//" from "//TRIM(filename)
 
   ! Some preparations
   obsfld=0.
@@ -152,20 +181,26 @@ program prep_obs
 
 
   ! User info - ADAPT ACCORDINGLY
-  write (*,*) " (prep_obs) Using data >40N and <45S"
+  WRITE (*,*) " (prep_obs) Using data >40N and <45S"
 
   ! Loop over space
-  do i=1,size(mask,1)
-    do j=1,size(mask,2)
+  DO i = 1, SIZE(mask, 1)
+    DO j = 1, SIZE(mask, 2)
 
       ! Pick out ocean points where data is existent
-      if ( (errorfld(i,j)) ** 2 .GT. 0.0001 .and. (errorfld(i,j)) ** 2 .LT. 1.0 &
-           .and. obsfld(i,j) .LT. 10.0 .and. obsfld(i,j) .GT. 0.0001 ) then
+      IF (       (errorfld(i,j)) ** 2 .GT. errmin ** 2  &
+           .AND. (errorfld(i,j)) ** 2 .LT. errmax ** 2  &
+           .AND. obsfld(i,j) .GT. obsmin                &
+           .AND. obsfld(i,j) .LT. obsmax ) THEN
         
-        ! Limit 'obs' to >40°N or <45°S
-        !if ( (lats(i,j).ge.40) .or. (lats(i,j).le.-45) ) then
-         if ( (lats(i,j) .ge. 40.0) ) then !.and. (lats(i,j) .le. 80.05) .and. (lons(i,j) .ge. 100.0)) then
-        !if ( (lats(i,j) .ge. 63.0)  .and. (lats(i,j) .le. 65.0))   then
+        ! Limit 'obs' spatially
+        IF (     (       lats(i,j) .GE. latmin_NH        & 
+                  .AND.  lats(i,j) .LE. latmax_NH      ) &
+             .OR.(       lats(i,j) .LE. (-latmin_SH)     &
+                  .AND.  lats(i,j) .GE. (-latmax_SH)   ) &
+           )  THEN 
+
+
           obscnt = obscnt + 1
           obs(obscnt)%d    = obsfld(i,j)
           obs(obscnt)%lon  = lons(i,j)
@@ -189,40 +224,35 @@ program prep_obs
           obs(obscnt)%date        = rdate
           obs(obscnt)%orig_id     =  0
 
-        end if ! >40°N or <50°S
-      end if ! ocean point
-    end do ! j
-  end do ! i
+        END IF ! Spatial selection
+      END IF ! if valid point
+    END DO ! j
+  END DO ! i
 
-  !obscnt = 3
   !Write data:
-  inquire(iolength=i)obs(1)
-  open (11, file='observations.uf', status='replace',form='unformatted', access='direct', recl=i)
-  do j=1,obscnt
-     write(11,rec=j)obs(j)
-  enddo
-  close(11)
-
-
-!  !XXX: 
-!  ! Write data to textfile, for visual checking
-  open(12, file='observations.txt')
-  do j=1,obscnt
-     write(12,FMT=53) obs(j)
+  INQUIRE(iolength=i)obs(1)
+  OPEN (11, file='observations.uf', status='replace',form='unformatted', access='direct', recl=i)
+  DO j = 1, obscnt
+     WRITE(11, rec=j)obs(j)
+  ENDDO
+  CLOSE(11)
+
+  ! Write data to textfile, for visual checking
+  OPEN(12, file = 'observations.txt')
+  DO j = 1, obscnt
+     WRITE(12, FMT = 53) obs(j)
 53   FORMAT(f8.4,X,f8.4,X,a8,X,2(f10.5,X),f4.2,X,2(I3,X),I1,X,4(f5.2,X),L,X,2(I3,X),f5.2,X,I8,X,I1)
-!      write(12, *) obs(j)
-  enddo
-  close(12)
-!  !:XXX
+  ENDDO
+  CLOSE(12)
 
   ! Tell user how many obs there were
-  write(*,*) " (prep_obs) Wrote out this many obs: ", obscnt
-  write(*,*) " (prep_obs) Number of ocean points : ", sum(mask)
+  WRITE(*,*) " (prep_obs) Wrote out this many obs: ", obscnt
+  WRITE(*,*) " (prep_obs) Number of ocean points : ", sum(mask)
 
   ! Cleanup
-  if (allocated(obs)) deallocate(obs,STAT=ierr)
+  IF (allocated(obs)) deallocate(obs,STAT=ierr)
 
-  write(*,*) ' (prep_obs) End successfully reached'; write(*,*)
+  WRITE(*,*) ' (prep_obs) End successfully reached'; WRITE(*,*)
 
 
 contains