prismtrs_apply_weights_2Dcons.F90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------
00002 ! Copyright 2007-2010, CERFACS, Toulouse, France.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PRISMTrs_Apply_weights_2Dcons
00008 !
00009 ! !INTERFACE
00010 subroutine prismtrs_apply_weights_2dcons(id_trans_out_size,   &
00011                                   dda_trans_out,       &
00012                                   id_trans_in_size,    &
00013                                   dda_trans_in,        &
00014                                   ida_mask,            &
00015                                   id_nb_links,         &  
00016                                   ida_grid1_add,       &
00017                                   ida_grid2_add,       &
00018                                   id_num_wts,          &
00019                                   dda_wts_map1,        &
00020                                   nbr_fields,          &
00021                                   id_err)
00022 !
00023 ! !USES:
00024 !
00025   USE PRISMDrv, dummy_interface => PRISMTrs_Apply_weights_2Dcons
00026   IMPLICIT NONE
00027 !
00028 ! !PARAMETERS:
00029 !
00030 ! size of the source transient
00031   INTEGER, INTENT (IN)                              :: id_trans_out_size
00032 
00033 ! Nbr of bundle components
00034   INTEGER, INTENT (IN)                              :: nbr_fields
00035 
00036 ! source transient field
00037   DOUBLE PRECISION, DIMENSION(id_trans_out_size*nbr_fields), INTENT(IN) :: dda_trans_out
00038  
00039 ! size of the target transient
00040   INTEGER, INTENT (IN)                              :: id_trans_in_size
00041 
00042 ! target transient mask
00043   INTEGER, DIMENSION(id_trans_in_size), INTENT(IN)       :: ida_mask
00044   
00045 ! number of links involved in the regridding
00046   INTEGER, INTENT (IN)                              :: id_nb_links
00047  
00048 ! source epio address for the links
00049   INTEGER, DIMENSION(id_nb_links), INTENT (IN)      :: ida_grid1_add
00050 
00051 ! target epio address for the links
00052   INTEGER, DIMENSION(id_nb_links), INTENT (IN)      :: ida_grid2_add
00053 
00054 ! Number of weights for 2D conservative remapping
00055   INTEGER, INTENT(in)                               :: id_num_wts
00056 
00057 ! weight for the links
00058   DOUBLE PRECISION, DIMENSION(id_num_wts,id_nb_links), INTENT(IN) :: dda_wts_map1
00059 
00060 !
00061 ! ! RETURN VALUE
00062 !
00063 ! transient out field
00064   DOUBLE PRECISION, DIMENSION(id_trans_in_size*nbr_fields), INTENT (Out) :: dda_trans_in
00065 
00066   INTEGER, INTENT (Out)               :: id_err   ! error value
00067 !
00068 ! !LOCAL VARIABLES
00069 !
00070 
00071 ! epio index of currently processed source cell
00072   INTEGER :: src_cell_idx
00073 ! epio index of currently/previous processed target cell
00074   INTEGER :: tgt_cell_idx, prev_tgt_cell_idx
00075 ! stride used for addressing each single field in the multiple field case
00076   INTEGER :: stride_in, stride_out
00077 ! !DESCRIPTION
00078 ! Subroutine "PRISMTrs_Interp_Apply_weights_2Dcons" performs the 
00079 ! matrix vector computation to get the field on the target grid 
00080 ! for 2D conservative remapping.
00081 !
00082 ! !REVISED HISTORY
00083 !   Date      Programmer   Description
00084 ! ----------  ----------   -----------
00085 ! 16/02/2007  S. Valcke    Creation
00086 ! 06/03/2008  J. Charles   Modifications added for the use of bundle fields
00087 !
00088 ! EOP
00089 !----------------------------------------------------------------------
00090 ! $Id: prismtrs_apply_weights_2Dcons.F90,v 1.1.2.7 2008/06/04 15:48:51 redler Exp $
00091 ! $Author: redler $
00092 !----------------------------------------------------------------------
00093 !
00094 ! Local declarations
00095 !
00096   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00097      '$Id: prismtrs_apply_weights_2Dcons.F90,v 1.1.2.7 2008/06/04 15:48:51 redler Exp $'
00098 !
00099   DOUBLE PRECISION :: 
00100      moy_tmp_S, moy_tmp_N,    ! field average at the pole
00101      moy_err_S, moy_err_N,    ! error average at the pole
00102      latpol_N, latpol_S,   
00103      dd_dif
00104 !
00105   DOUBLE PRECISION :: dst_err(id_trans_in_size)
00106 ! 
00107   INTEGER :: ib, ib_bis, i
00108   INTEGER :: compt_S, compt_N  ! number of cells representating a pole
00109 !
00110 ! ---------------------------------------------------------------------
00111 !
00112 #ifdef VERBOSE
00113   PRINT *, '| | | | | | | Enter PRISMTrs_Apply_weights_2Dcons'
00114   call psmile_flushstd
00115 #endif
00116 !
00117   id_err = 0
00118 !
00119 ! 1. For each non masked target point, apply the weights on the source field 
00120 !    to get the target field
00121 !
00122    dda_trans_in(:) = psmile_dundef
00123    stride_in = 0
00124    stride_out = 0
00125    tgt_cell_idx = -1
00126 
00127    DO i = 1, nbr_fields
00128       DO ib_bis = 1, id_nb_links, 5000
00129          DO ib = ib_bis, min(id_nb_links, ib_bis+5000-1)
00130             prev_tgt_cell_idx = tgt_cell_idx
00131             ! get the epio indices for the source and target cell of the current source-target-cell-link
00132             tgt_cell_idx = ida_grid2_add(ib) + stride_in
00133             src_cell_idx = ida_grid1_add(ib) + stride_out
00134 #ifdef PRISM_ASSERTION
00135             ! the target epio address array needs to be sorted. this is required for detecting target
00136             ! cell which have no corresponding source cell
00137             ! the sorting is done in prismtrs_prismtrs_conserv2d_weight.F90 using the routine
00138             ! prismtrs_sort_add
00139             if (tgt_cell_idx < prev_tgt_cell_idx) then
00140                print *, trim(ch_id), 'tgt_cell_idx, prev_tgt_cell_idx', &
00141                                        tgt_cell_idx, prev_tgt_cell_idx
00142                ! TODO: psmile_assert should be moved to psmile_common in order be usable within the driver/transformer
00143 !              call psmile_assert (__FILE__, __LINE__, &
00144 !                                 "source epio address array is not sorted correctly")
00145                print *, "Assertion violation: source epio address array is not sorted correctly"
00146                print *, "file", __FILE__, "line", __LINE__ - 7
00147                call psmile_flushstd
00148                call psmile_abort
00149             endif
00150 #endif
00151 #ifdef DEBUGX
00152    IF(ida_grid2_add(ib) == 7447) WRITE(ida_grid2_add(ib),*) ida_grid2_add(ib), ida_grid1_add(ib),dda_wts_map1(1,ib)
00153 #endif
00154             ! if the current target cell is touched for the first time => init with zero
00155             IF (prev_tgt_cell_idx .ne. tgt_cell_idx) dda_trans_in(tgt_cell_idx) = 0.0
00156             ! apply the weight for the current source-target-cell-link
00157             dda_trans_in(tgt_cell_idx) = dda_trans_in(tgt_cell_idx) &
00158                   + dda_wts_map1(1,ib) * dda_trans_out(src_cell_idx)
00159          ENDDO
00160       ENDDO
00161       ! set the strides for the next field
00162       stride_in  = stride_in  + id_trans_in_size
00163       stride_out = stride_out + id_trans_out_size
00164    ENDDO
00165   !
00166 !!$  !    Check for line of points at the poles (as done in scriprmp.F)
00167 !!$  !    ------------------------------------------------------------
00168 !!$  !
00169 !!$  !    For the north pole _N and the south pole _S
00170 !!$
00171 !!$  latpol_N = dble_pi*half
00172 !!$  latpol_S = -dble_pi*half
00173 !!$  compt_N = 0
00174 !!$  moy_tmp_N = 0d0
00175 !!$  moy_err_N = 0d0
00176 !!$  compt_S = 0
00177 !!$  moy_tmp_S = 0d0
00178 !!$  moy_err_S = 0d0
00179 !!$  !
00180 !!$  DO ib = 1, id_trans_in_size
00181 !!$    IF (dst_lat(ib) == latpol_N) THEN
00182 !!$        moy_tmp_N = moy_tmp_N + dda_trans_in(ib)
00183 !!$        compt_N = compt_N + 1
00184 !!$    ELSE IF (dst_lat(ib) == latpol_S) THEN
00185 !!$        moy_tmp_S = moy_tmp_S + dda_trans_in(ib)
00186 !!$        compt_S = compt_S + 1
00187 !!$    END IF
00188 !!$  ENDDO
00189 !!$  !
00190 !!$  IF (compt_N/=0) THEN
00191 !!$      moy_tmp_N = moy_tmp_N/compt_N
00192 !!$      DO ib = 1, id_trans_in_size
00193 !!$        IF (dst_lat(ib) == latpol_N) THEN
00194 !!$            dst_err(ib) = abs(moy_tmp_N - dda_trans_in(ib))
00195 !!$            moy_err_N = moy_err_N + dst_err(ib)
00196 !!$            dda_trans_in(ib) = moy_tmp_N
00197 !!$        END IF
00198 !!$      END DO
00199 !!$      moy_err_N = moy_err_N/compt_N
00200 !!$      PRINT*, 'Points at the pole N'
00201 !!$      PRINT*, 'Average value at the pole : ', moy_tmp_N
00202 !!$      PRINT*, 'Average error at the pole : ', moy_err_N
00203 !!$  ELSE
00204 !!$      PRINT*, 'No point at the pole N'
00205 !!$  END IF
00206 !!$  !
00207 !!$  IF (compt_S/=0) THEN
00208 !!$      moy_tmp_S = moy_tmp_S/compt_S
00209 !!$      DO ib = 1, id_trans_in_size
00210 !!$        IF (dst_lat(ib) == latpol_S) THEN
00211 !!$            dst_err(ib) = ABS(moy_tmp_S - dda_trans_in(ib))
00212 !!$            moy_err_S = moy_err_S + dst_err(ib)
00213 !!$            dda_trans_in(ib) = moy_tmp_S
00214 !!$        END IF
00215 !!$      END DO
00216 !!$      moy_err_S = moy_err_S/compt_S
00217 !!$      PRINT*, 'Points at the pole S'
00218 !!$      PRINT*, 'Average value at the pole : ', moy_tmp_S
00219 !!$      PRINT*, 'Average error at the pole : ', moy_err_S
00220 !!$  ELSE
00221 !!$      PRINT*, 'No point at the pole S'
00222 !!$  END IF
00223 !
00224 #ifdef VERBOSE
00225   PRINT *, '| | | | | | | Quit PRISMTrs_Apply_weights_2Dcons'
00226   call psmile_flushstd
00227 #endif
00228 
00229 END SUBROUTINE PRISMTrs_Apply_weights_2Dcons

Generated on 18 Mar 2011 for Oasis4 by  doxygen 1.6.1