Oasis3 4.0.2
|
00001 SUBROUTINE grstat (pcorx, pcory, psgr, kmsk, kngx, kngy, cdtyp, 00002 $ pdr, psr, phhi, kvma) 00003 C**** 00004 C ***************************** 00005 C * OASIS ROUTINE - LEVEL T * 00006 C * ------------- ------- * 00007 C ***************************** 00008 C 00009 C**** *grstat* - Statistic routine 00010 C 00011 C Purpose: 00012 C ------- 00013 C Given a grid, calculate the average grid square size, 00014 C the total surface and the inverse of sum of surface elements squared. 00015 C 00016 C** Interface: 00017 C --------- 00018 C *CALL* *grstat (pcorx, pcory, psgr, kmsk, kngx, kngy, cdtyp, 00019 C pdr, psr, phhi, kvma)* 00020 C 00021 C Input: 00022 C ----- 00023 C pcorx : the grid points longitude (real 2D) 00024 C pcory : the grid points latitude (real 2D) 00025 C psgr : surface elements (real 2D) 00026 C kmsk : the grid mask (integer 2D) 00027 C kngx : the grid size in direction 1 00028 C kngy : the grid size in direction 2 00029 C cdtyp : grid type 00030 C cdtyp : type of source grid 00031 C kvma : integer value of the mask 00032 C 00033 C Output: 00034 C ------ 00035 C pdr : the average grid square size 00036 C psr : the total surface 00037 C phhi : inverse of sum of surface elements squared 00038 C 00039 C Workspace: 00040 C --------- 00041 C None 00042 C 00043 C Externals: 00044 C --------- 00045 C sqdis 00046 C 00047 C Reference: 00048 C --------- 00049 C O. Thual, Simple ocean-atmosphere interpolation. 00050 C Part A: The method, EPICOA 0629 (1992) 00051 C Part B: Software implementation, EPICOA 0630 (1992) 00052 C 00053 C See also OASIS manual (1995) 00054 C 00055 C History: 00056 C ------- 00057 C Version Programmer Date Description 00058 C ------- ---------- ---- ----------- 00059 C 1.1 O. Thual 94/01/01 created 00060 C 2.0 L. Terray 95/12/26 modified : to suit OASIS 2.0 00061 C 00062 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00063 C 00064 C* ---------------------------- Include files --------------------------- 00065 C 00066 USE mod_printing 00067 USE mod_unit 00068 C 00069 C* ---------------------------- Argument declarations ------------------- 00070 C 00071 REAL (kind=ip_realwp_p) pcorx(kngx,kngy), pcory(kngx,kngy) 00072 REAL (kind=ip_realwp_p) psgr(kngx,kngy) 00073 INTEGER (kind=ip_intwp_p) kmsk(kngx,kngy) 00074 CHARACTER*1 cdtyp 00075 C 00076 C* ---------------------------- Local declarations ---------------------- 00077 C 00078 CHARACTER*1 cltyp 00079 C 00080 C* ---------------------------- Poema verses ---------------------------- 00081 C 00082 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00083 C 00084 C* 1. Initialization 00085 C -------------- 00086 C 00087 IF (nlogprt .GE. 2) THEN 00088 WRITE (UNIT = nulou,FMT = *) ' ' 00089 WRITE (UNIT = nulou,FMT = *) '***** * Entering ROUTINE grstat' 00090 WRITE (UNIT = nulou,FMT = *) '*******************************' 00091 WRITE (UNIT = nulou,FMT = *) ' ' 00092 call flush(nulou) 00093 ENDIF 00094 C 00095 ivma = kvma 00096 cltyp = cdtyp 00097 C 00098 C 00099 C* 2. Interdistances only for unmasked points 00100 C --------------------------------------- 00101 C 00102 C* Average distance squared in x-direction 00103 C 00104 ztot2 = 0. 00105 icount = 0 00106 DO 210 j2 = 1, kngy 00107 zsum = 0. 00108 DO 220 j1 = 1, kngx-1 00109 IF (kmsk(j1,j2) .NE. ivma) THEN 00110 iadd = 0 00111 IF (cltyp .ne. 'D') iadd = 1 00112 C* For reduced grid, calculate distance only in the x-direction between 00113 C* points located on same latitude circle. 00114 IF (cltyp .eq. 'D' .and. pcory(j1,j2) .eq. pcory(j1+1,j2)) 00115 $ iadd = 1 00116 C* For U grid, the average distance will be calculated based on all 00117 C* consecutive points in the arrays, which is not exact 00118 IF (iadd .eq. 1) then 00119 zdis = sqdis(pcorx(j1,j2), pcory(j1,j2), 00120 $ pcorx(j1+1,j2), pcory(j1+1,j2)) 00121 icount = icount + iadd 00122 zsum = zsum + zdis 00123 ENDIF 00124 ENDIF 00125 220 CONTINUE 00126 ztot2 = ztot2 + zsum 00127 210 CONTINUE 00128 WRITE(nulan, *)'In grstat, icount= ', icount 00129 ztot2 = ztot2 / float(icount) 00130 C 00131 C* Average distance squared in y-direction 00132 C 00133 IF ((cltyp .eq. 'D') .or. (cltyp .eq. 'U')) THEN 00134 pdr = ztot2 00135 ztot1 = 0. 00136 ELSE 00137 ztot1 = 0. 00138 icount = 0 00139 DO 230 j1 = 1, kngx 00140 zsum = 0. 00141 DO 240 j2 = 1, kngy-1 00142 IF(kmsk(j1,j2) .NE. ivma) THEN 00143 zdis = sqdis(pcorx(j1,j2), pcory(j1,j2), 00144 $ pcorx(j1,j2+1), pcory(j1,j2+1)) 00145 icount = icount + 1 00146 zsum = zsum + zdis 00147 ENDIF 00148 240 CONTINUE 00149 ztot1 = ztot1 + zsum 00150 230 CONTINUE 00151 ztot1 = ztot1 / float(icount) 00152 C 00153 C* Get average grid square size 00154 C 00155 pdr = (ztot1 + ztot2) / 2. 00156 ENDIF 00157 00158 C 00159 C 00160 C* 3. Surface and <hh>^-1 00161 C ------------------- 00162 C 00163 zsurf = 0. 00164 zhhi = 0. 00165 DO 310 j2 = 1, kngy 00166 DO 320 j1 = 1, kngx 00167 IF (kmsk(j1,j2) .NE. ivma) THEN 00168 zsurf = zsurf + psgr(j1,j2) 00169 zhhi = zhhi + psgr(j1,j2) * psgr(j1,j2) 00170 ENDIF 00171 320 CONTINUE 00172 310 CONTINUE 00173 C 00174 C* Get total surface and inverse of sum of elements squared 00175 C 00176 psr = zsurf 00177 IF (zhhi .ne. 0.) phhi = 1. / zhhi 00178 C 00179 C 00180 C* 4. End of routine 00181 C -------------- 00182 C 00183 IF (nlogprt .GE. 2) THEN 00184 WRITE (UNIT = nulou,FMT = *) ' ' 00185 WRITE (UNIT = nulou,FMT = *) '******** Leaving ROUTINE grstat' 00186 WRITE (UNIT = nulou,FMT = *) '*******************************' 00187 WRITE (UNIT = nulou,FMT = *) ' ' 00188 call flush(nulou) 00189 ENDIF 00190 RETURN 00191 END