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 00135 00136 00137 #ifdef PRISM_ASSERTION 00138 ! the target epio address array needs to be sorted. this is required for detecting target 00139 ! cell which have no corresponding source cell 00140 ! the sorting is done in prismtrs_prismtrs_conserv2d_weight.F90 using the routine 00141 ! prismtrs_sort_add 00142 if (tgt_cell_idx < prev_tgt_cell_idx) then 00143 print *, trim(ch_id), 'tgt_cell_idx, prev_tgt_cell_idx', & 00144 tgt_cell_idx, prev_tgt_cell_idx 00145 ! TODO: psmile_assert should be moved to psmile_common in order be usable within the driver/transformer 00146 ! call psmile_assert (__FILE__, __LINE__, & 00147 ! "source epio address array is not sorted correctly") 00148 print *, "Assertion violation: source epio address array is not sorted correctly" 00149 print *, "file", __FILE__, "line", __LINE__ - 7 00150 call psmile_flushstd 00151 call psmile_abort 00152 endif 00153 #endif 00154 00155 #ifdef DEBUGX 00156 IF(ida_grid2_add(ib) == 2465) WRITE(88,*) 'In apply_weights_2Dcons' 00157 IF(ida_grid2_add(ib) == 2465) WRITE(88,*) ib,ida_grid2_add(ib), ida_grid1_add(ib),dda_wts_map1(1,ib) 00158 call psmile_flushstd 00159 #endif 00160 ! if the current target cell is touched for the first time => init with zero 00161 IF (prev_tgt_cell_idx .ne. tgt_cell_idx) dda_trans_in(tgt_cell_idx) = 0.0 00162 ! apply the weight for the current source-target-cell-link 00163 dda_trans_in(tgt_cell_idx) = dda_trans_in(tgt_cell_idx) & 00164 + dda_wts_map1(1,ib) * dda_trans_out(src_cell_idx) 00165 ENDDO 00166 ENDDO 00167 ! set the strides for the next field 00168 stride_in = stride_in + id_trans_in_size 00169 stride_out = stride_out + id_trans_out_size 00170 ENDDO 00171 ! 00172 !!$ ! Check for line of points at the poles (as done in scriprmp.F) 00173 !!$ ! ------------------------------------------------------------ 00174 !!$ ! 00175 !!$ ! For the north pole _N and the south pole _S 00176 !!$ 00177 !!$ latpol_N = dble_pi*half 00178 !!$ latpol_S = -dble_pi*half 00179 !!$ compt_N = 0 00180 !!$ moy_tmp_N = 0d0 00181 !!$ moy_err_N = 0d0 00182 !!$ compt_S = 0 00183 !!$ moy_tmp_S = 0d0 00184 !!$ moy_err_S = 0d0 00185 !!$ ! 00186 !!$ DO ib = 1, id_trans_in_size 00187 !!$ IF (dst_lat(ib) == latpol_N) THEN 00188 !!$ moy_tmp_N = moy_tmp_N + dda_trans_in(ib) 00189 !!$ compt_N = compt_N + 1 00190 !!$ ELSE IF (dst_lat(ib) == latpol_S) THEN 00191 !!$ moy_tmp_S = moy_tmp_S + dda_trans_in(ib) 00192 !!$ compt_S = compt_S + 1 00193 !!$ END IF 00194 !!$ ENDDO 00195 !!$ ! 00196 !!$ IF (compt_N/=0) THEN 00197 !!$ moy_tmp_N = moy_tmp_N/compt_N 00198 !!$ DO ib = 1, id_trans_in_size 00199 !!$ IF (dst_lat(ib) == latpol_N) THEN 00200 !!$ dst_err(ib) = abs(moy_tmp_N - dda_trans_in(ib)) 00201 !!$ moy_err_N = moy_err_N + dst_err(ib) 00202 !!$ dda_trans_in(ib) = moy_tmp_N 00203 !!$ END IF 00204 !!$ END DO 00205 !!$ moy_err_N = moy_err_N/compt_N 00206 !!$ PRINT*, 'Points at the pole N' 00207 !!$ PRINT*, 'Average value at the pole : ', moy_tmp_N 00208 !!$ PRINT*, 'Average error at the pole : ', moy_err_N 00209 !!$ ELSE 00210 !!$ PRINT*, 'No point at the pole N' 00211 !!$ END IF 00212 !!$ ! 00213 !!$ IF (compt_S/=0) THEN 00214 !!$ moy_tmp_S = moy_tmp_S/compt_S 00215 !!$ DO ib = 1, id_trans_in_size 00216 !!$ IF (dst_lat(ib) == latpol_S) THEN 00217 !!$ dst_err(ib) = ABS(moy_tmp_S - dda_trans_in(ib)) 00218 !!$ moy_err_S = moy_err_S + dst_err(ib) 00219 !!$ dda_trans_in(ib) = moy_tmp_S 00220 !!$ END IF 00221 !!$ END DO 00222 !!$ moy_err_S = moy_err_S/compt_S 00223 !!$ PRINT*, 'Points at the pole S' 00224 !!$ PRINT*, 'Average value at the pole : ', moy_tmp_S 00225 !!$ PRINT*, 'Average error at the pole : ', moy_err_S 00226 !!$ ELSE 00227 !!$ PRINT*, 'No point at the pole S' 00228 !!$ END IF 00229 ! 00230 #ifdef VERBOSE 00231 PRINT *, '| | | | | | | Quit PRISMTrs_Apply_weights_2Dcons' 00232 call psmile_flushstd 00233 #endif 00234 00235 END SUBROUTINE PRISMTrs_Apply_weights_2Dcons
1.6.1