Oasis3 4.0.2
|
00001 SUBROUTINE pcssph (px, py, psgr, kngx, kngy) 00002 C**** 00003 C ***************************** 00004 C * OASIS ROUTINE - LEVEL T * 00005 C * ------------- ------- * 00006 C ***************************** 00007 C 00008 C**** *pcssph* - Arithmetic routine 00009 C 00010 C Purpose: 00011 C ------- 00012 C Calculate surface element for a spheric and periodic grid. 00013 C The coordinates are in degrees, there are no pole points. 00014 C 00015 C** Interface: 00016 C --------- 00017 C *CALL* *pcssph(px, py, psgr, kngx, kngy)* 00018 C 00019 C Input: 00020 C ----- 00021 C px : grid longitudes (real 2D) 00022 C py : grid latitudes (real 2D) 00023 C kngx : number of longitudes 00024 C kngy : number of latitudes 00025 C 00026 C Output: 00027 C ------ 00028 C psgr : grid surface elements (real 2D) 00029 C 00030 C 00031 C Workspace: 00032 C --------- 00033 C None 00034 C 00035 C External: 00036 C -------- 00037 C None 00038 C 00039 C References: 00040 C ---------- 00041 C O. Thual, Simple ocean-atmosphere interpolation. 00042 C Part A: The method, EPICOA 0629 (1992) 00043 C Part B: Software implementation, EPICOA 0630 (1992) 00044 C See also OASIS manual (1995) 00045 C 00046 C History: 00047 C ------- 00048 C Version Programmer Date Description 00049 C ------- ---------- ---- ----------- 00050 C 1.1 O. Thual 93/04/15 created 00051 C 2.0 L. Terray 95/10/01 modified: new structure 00052 C 00053 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00054 C 00055 C* ---------------------------- Include files --------------------------- 00056 C 00057 USE mod_unit 00058 USE mod_printing 00059 C 00060 C* ---------------------------- Argument declarations ------------------- 00061 C 00062 REAL (kind=ip_realwp_p) px(kngx,kngy), py(kngx,kngy) 00063 REAL (kind=ip_realwp_p) psgr(kngx,kngy) 00064 C 00065 C* ---------------------------- Poema verses ---------------------------- 00066 C 00067 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00068 C 00069 IF (nlogprt .GE. 2) THEN 00070 WRITE(nulou,*) '***** * Entering ROUTINE pcssph ******* ' 00071 CALL FLUSH(nulou) 00072 ENDIF 00073 C 00074 C* 1. Degres to radians conversion factor 00075 C ----------------------------------- 00076 C 00077 zconv = 1.74532925199432957692e-2 00078 C 00079 C 00080 C* 2. Surfaces 00081 C -------- 00082 C 00083 DO 210 j2 = 1, kngy 00084 DO 220 j1 = 1, kngx 00085 C 00086 C* Left or right periodicity 00087 C 00088 IF (J1 .EQ. 1) THEN 00089 zx1 = px(kngx,j2) - 360. 00090 zx2 = px(2,j2) 00091 ELSE IF (j1 .EQ. kngx) THEN 00092 zx1 = px(kngx-1,j2) 00093 zx2 = px(1,j2) + 360. 00094 ELSE 00095 zx1 = px(j1-1,j2) 00096 zx2 = px(j1+1,j2) 00097 ENDIF 00098 C 00099 C* Bottom or top treatment 00100 C 00101 IF (j2 .EQ. 1) THEN 00102 zy1 = -90. 00103 zy2 = .5 * (py(j1,j2) + py(j1,j2+1)) 00104 ELSE IF (j2 .EQ. kngy) THEN 00105 zy1 = .5 * (py(j1,j2) + py(j1,j2-1)) 00106 zy2 = 90. 00107 ELSE 00108 zy1 = .5 * (py(j1,j2) + py(j1,j2-1)) 00109 zy2 = .5 * (py(j1,j2) + py(j1,j2+1)) 00110 ENDIF 00111 C 00112 C* Conversion to radians 00113 C 00114 zfi1 = zx1 * zconv 00115 zfi2 = zx2 * zconv 00116 zth1 = zy1 * zconv 00117 zth2 = zy2 * zconv 00118 C 00119 C* Calculate grid square surface 00120 C 00121 zfac = sin(zth2) - sin(zth1) 00122 psgr(j1,j2) = abs(.5 * (zfi2-zfi1) * zfac) 00123 220 continue 00124 210 CONTINUE 00125 C 00126 IF (nlogprt .GE. 2) THEN 00127 WRITE(nulou,*) '***** * Leaving ROUTINE pcssph ******* ' 00128 CALL FLUSH(nulou) 00129 ENDIF 00130 C* End of routine 00131 C 00132 RETURN 00133 END 00134 00135 00136 00137 00138 00139 00140