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