|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE SSO(OSSO,OSSO_ANIS,PSEA) 00003 ! ######################### 00004 ! 00005 !!*SSO computes the SSO anisotropy, direction and slope 00006 !! 00007 !! 00008 !! METHOD 00009 !! ------ 00010 !! See Lott and Miller, 1997, QJRMS 101-127 00011 !! 00012 !! AUTHOR 00013 !! ------ 00014 !! 00015 !! V.Masson Meteo-France 00016 !! 00017 !! MODIFICATION 00018 !! ------------ 00019 !! 00020 !! Original 23/05/97 00021 !! 00022 !---------------------------------------------------------------------------- 00023 ! 00024 !* 0. DECLARATION 00025 ! ----------- 00026 ! 00027 USE MODI_GET_MESH_DIM 00028 USE MODI_GET_ADJACENT_MESHES 00029 ! 00030 USE MODD_CSTS, ONLY : XPI 00031 USE MODD_PGDWORK, ONLY : NSSO, XSSQO, LSSQO 00032 USE MODD_PGD_GRID, ONLY : NL, CGRID, XGRID_PAR, NGRID_PAR 00033 USE MODD_SURF_ATM_GRID_n, ONLY : XMESH_SIZE 00034 USE MODD_SURF_ATM_SSO_n, ONLY : XSSO_ANIS, XSSO_DIR, XSSO_SLOPE 00035 ! 00036 ! 00037 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00038 USE PARKIND1 ,ONLY : JPRB 00039 ! 00040 IMPLICIT NONE 00041 ! 00042 !* 0.1 Declaration of dummy arguments 00043 ! ------------------------------ 00044 ! 00045 LOGICAL, DIMENSION(:), INTENT(OUT) :: OSSO ! .T. : the SSO coefficients 00046 ! ! are computed at grid point 00047 ! ! .F. : not enough sub-grid 00048 ! ! information avalaible to 00049 ! ! compute the coefficients 00050 LOGICAL, DIMENSION(:), INTENT(OUT) :: OSSO_ANIS ! .T. : the SSO anisotropy 00051 ! ! are computed at grid point 00052 ! ! .F. : not enough sub-grid 00053 ! ! information avalaible to 00054 ! ! compute the coefficients 00055 REAL, DIMENSION(:), INTENT(IN) :: PSEA ! sea fraction 00056 ! 00057 !* 0.2 Declaration of indexes 00058 ! ---------------------- 00059 ! 00060 INTEGER :: JL ! loop index on grid meshs 00061 INTEGER :: IL ! grid mesh index of second subgrid point used 00062 INTEGER :: JISS, JJSS ! loop indexes for subsquares arrays 00063 INTEGER :: JNEXT ! loop index on subgrid meshes 00064 INTEGER :: JPREV ! loop index on subgrid meshes 00065 INTEGER :: INEXT ! index to add to JISS or JJSS to obtain the following 00066 ! ! point containing a data in a segment 00067 INTEGER :: IPREV ! index to remove from JISS or JJSS to obtain the 00068 ! ! previous point containing a data in a segment 00069 ! 00070 INTEGER :: IMAXI ! index of the next subsquare with data along I axis, 00071 ! or last subsquare inside grid mesh along I axis 00072 INTEGER :: IMAXJ ! index of the next subsquare with data along J axis, 00073 ! or last subsquare inside grid mesh along J axis 00074 INTEGER, DIMENSION(NL) :: ILEFT ! index of left grid mesh 00075 INTEGER, DIMENSION(NL) :: IRIGHT ! index of right grid mesh 00076 INTEGER, DIMENSION(NL) :: ITOP ! index of top grid mesh 00077 INTEGER, DIMENSION(NL) :: IBOTTOM ! index of bottom grid mesh 00078 ! 00079 !* 0.3 Declaration of working arrays inside a MESONH grid (JI,JJ) 00080 ! ----------------------------- 00081 ! 00082 REAL, DIMENSION(NSSO,NSSO) :: ZDZSDX ! topographic gradient along x 00083 REAL, DIMENSION(NSSO,NSSO) :: ZDZSDY ! topographic gradient along y 00084 LOGICAL, DIMENSION(NSSO,NSSO) :: GDZSDX ! occurence of gradient along x 00085 LOGICAL, DIMENSION(NSSO,NSSO) :: GDZSDY ! occurence of gradient along y 00086 REAL :: ZDXEFF ! width of a subsquare along I axis 00087 REAL :: ZDYEFF ! width of a subsquare along J axis 00088 LOGICAL :: GBOUND ! .T.: first left (for dzs/dx) or first bottom (for dzs/dy) 00089 ! sub-square gradients is being computed 00090 ! 00091 !* 0.4 Declaration of other local variables 00092 ! ------------------------------------ 00093 ! 00094 REAL, DIMENSION(NL) :: ZDX ! grid mesh size in x direction 00095 REAL, DIMENSION(NL) :: ZDY ! grid mesh size in y direction 00096 REAL, DIMENSION(NL) :: ZHXX ! topographic gradient correlation tensor: x,x 00097 REAL, DIMENSION(NL) :: ZHXY ! topographic gradient correlation tensor: x,y 00098 REAL, DIMENSION(NL) :: ZHYY ! topographic gradient correlation tensor: y,y 00099 REAL, DIMENSION(NL) :: ZK, ZL, ZM ! diagonalised terms of the tensor 00100 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00101 ! 00102 !---------------------------------------------------------------------------- 00103 ! 00104 !* 1. Initializations 00105 ! --------------- 00106 ! 00107 IF (LHOOK) CALL DR_HOOK('SSO',0,ZHOOK_HANDLE) 00108 ZK=0. 00109 ZL=0. 00110 ZM=0. 00111 ! 00112 !* 1.1 Occurence of computation of the coefficients 00113 ! -------------------------------------------- 00114 ! 00115 OSSO (:)=.FALSE. 00116 OSSO_ANIS(:)=.FALSE. 00117 ! 00118 ! 00119 !* 1.2 Grid dimension (meters) 00120 ! ----------------------- 00121 ! 00122 CALL GET_MESH_DIM(CGRID,NGRID_PAR,NL,XGRID_PAR,ZDX,ZDY,XMESH_SIZE) 00123 ! 00124 ! 00125 !* 1.3 Left, top, right and bottom adjacent gris meshes 00126 ! ------------------------------------------------ 00127 ! 00128 CALL GET_ADJACENT_MESHES(CGRID,NGRID_PAR,NL,XGRID_PAR,ILEFT,IRIGHT,ITOP,IBOTTOM) 00129 ! 00130 !---------------------------------------------------------------------------- 00131 ! 00132 !* 2. Loop on MESONH grid points 00133 ! -------------------------- 00134 ! 00135 DO JL=1,NL 00136 ! 00137 ! 00138 !* 2.1 No land in grid mesh 00139 ! -------------------- 00140 ! 00141 IF (PSEA(JL)==1.) CYCLE 00142 ! 00143 ZDXEFF=ZDX(JL)/FLOAT(NSSO) 00144 ZDYEFF=ZDY(JL)/FLOAT(NSSO) 00145 ! 00146 !* 2.3 Not enough data for computation 00147 ! ------------------------------- 00148 ! 00149 ! 1st step: removes points where orography data is not present at all. 00150 ! 00151 IF ( COUNT( LSSQO(:,:,JL)) == 0 ) CYCLE 00152 ! 00153 !---------------------------------------------------------------------------- 00154 ! 00155 !* 3. Computations of the gradients along x 00156 ! ------------------------------------- 00157 ! 00158 GDZSDX (:,:)=.FALSE. 00159 ZDZSDX (:,:)= 0. 00160 ! 00161 !* 3.1 loop on jss index (there is no specific computation along j) 00162 ! ----------------- 00163 ! 00164 DO JJSS=1,NSSO 00165 ! 00166 !* left point mark initialization 00167 ! 00168 GBOUND=.TRUE. 00169 ! 00170 !* 3.2 loop on iss index 00171 ! ----------------- 00172 ! 00173 DO JISS=1,NSSO 00174 ! 00175 ! 00176 !* 3.3 search for two consecutive grid points 00177 ! -------------------------------------- 00178 ! 00179 !* 3.3.1 first one 00180 ! 00181 IF (.NOT. LSSQO(JISS,JJSS,JL) ) CYCLE 00182 ! 00183 !* 3.3.2 second one (up to one grid mesh further) 00184 ! 00185 DO JNEXT=1,NSSO 00186 IF (JISS+JNEXT>NSSO) THEN 00187 IL = IRIGHT(JL) 00188 INEXT = JISS+JNEXT-NSSO 00189 ELSE 00190 IL = JL 00191 INEXT = JISS+JNEXT 00192 END IF 00193 ! no right point 00194 IF (IL==0) EXIT 00195 ! subgrid data found 00196 IF (LSSQO(INEXT,JJSS,IL)) EXIT 00197 END DO 00198 ! 00199 !* 3.3.3 none found: end of loop along jss 00200 ! --------------------- 00201 ! 00202 IF (JNEXT>=NSSO+1) EXIT 00203 ! 00204 !* 3.3.4 second point outside of the domain: end of loop along iss 00205 ! --------------------- 00206 ! 00207 IF (IL==0) EXIT 00208 ! 00209 !* 3.4 dzs/dx term 00210 ! ----------- 00211 ! 00212 IMAXI = MIN(JISS+JNEXT-1,NSSO) 00213 ! 00214 ZDZSDX(JISS:IMAXI,JJSS) = ( XSSQO(INEXT,JJSS,IL) - XSSQO(JISS,JJSS,JL)) & 00215 / FLOAT(JNEXT) / ZDXEFF 00216 ! 00217 GDZSDX(JISS:IMAXI,JJSS) = .TRUE. 00218 ! 00219 ! 00220 !* 3.5 left data point not on the left of the grid mesh (one more computation) 00221 ! ------------------------------------------------ 00222 ! 00223 IF (GBOUND .AND. JISS/=1) THEN 00224 DO JPREV=1,NSSO 00225 IF (JISS-JPREV<1) THEN 00226 IL = ILEFT(JL) 00227 IPREV = JISS-JPREV+NSSO 00228 ELSE 00229 IL = JL 00230 IPREV = JISS-JPREV 00231 END IF 00232 ! no left point 00233 IF (IL==0) EXIT 00234 ! subgrid data found 00235 IF (LSSQO(IPREV,JJSS,IL)) EXIT 00236 END DO 00237 00238 IF (.NOT. (JPREV>=NSSO+1 .OR. IL==0)) THEN 00239 ZDZSDX(1:JISS,JJSS) = ( XSSQO(JISS,JJSS,JL) - XSSQO(IPREV,JJSS,IL)) & 00240 / FLOAT(JPREV) / ZDXEFF 00241 ! 00242 GDZSDX(1:JISS,JJSS) = .TRUE. 00243 END IF 00244 END IF 00245 ! 00246 ! 00247 GBOUND=.FALSE. 00248 ! 00249 ! 00250 !* 3.6 end of loop on iss index 00251 ! ------------------------ 00252 ! 00253 END DO 00254 ! 00255 !* 3.7 end of loop on jss index 00256 ! ------------------------ 00257 ! 00258 END DO 00259 ! 00260 !---------------------------------------------------------------------------- 00261 ! 00262 !* 4. Computations of the gradients along y 00263 ! ------------------------------------- 00264 ! 00265 GDZSDY (:,:)=.FALSE. 00266 ZDZSDY (:,:)= 0. 00267 ! 00268 !* 4.1 loop on iss index (there is no specific computation along i) 00269 ! ----------------- 00270 ! 00271 DO JISS=1,NSSO 00272 ! 00273 !* bottom point mark initialization 00274 ! 00275 GBOUND=.TRUE. 00276 ! 00277 !* 4.2 loop on jss index 00278 ! ----------------- 00279 ! 00280 DO JJSS=1,NSSO 00281 ! 00282 ! 00283 !* 4.3 search for two consecutive grid points 00284 ! -------------------------------------- 00285 ! 00286 !* 4.3.1 first one 00287 ! 00288 IF (.NOT. LSSQO(JISS,JJSS,JL) ) CYCLE 00289 ! 00290 !* 4.3.2 second one (up to one grid mesh further) 00291 ! 00292 DO JNEXT=1,NSSO 00293 IF (JJSS+JNEXT>NSSO) THEN 00294 IL = ITOP(JL) 00295 INEXT = JJSS+JNEXT-NSSO 00296 ELSE 00297 IL = JL 00298 INEXT = JJSS+JNEXT 00299 END IF 00300 ! no top point 00301 IF (IL==0) EXIT 00302 ! subgrid data found 00303 IF (LSSQO(JISS,INEXT,IL)) EXIT 00304 END DO 00305 ! 00306 !* 4.3.3 none found: end of loop along iss 00307 ! --------------------- 00308 ! 00309 IF (JNEXT>=NSSO+1) EXIT 00310 ! 00311 !* 3.3.4 second point outside of the domain: end of loop along iss 00312 ! --------------------- 00313 ! 00314 IF (IL==0) EXIT 00315 ! 00316 !* 4.4 dzs/dy term 00317 ! ----------- 00318 ! 00319 IMAXJ = MIN(JJSS+JNEXT-1,NSSO) 00320 ! 00321 ZDZSDY(JISS,JJSS:IMAXJ) = ( XSSQO(JISS,INEXT,IL) - XSSQO(JISS,JJSS,JL)) & 00322 / FLOAT(JNEXT) / ZDYEFF 00323 ! 00324 GDZSDY(JISS,JJSS:IMAXJ) = .TRUE. 00325 ! 00326 ! 00327 !* 4.5 bottom data point not on the bottom of the grid mesh (one more computation) 00328 ! ---------------------------------------------------- 00329 ! 00330 IF (GBOUND .AND. JJSS/=1) THEN 00331 DO JPREV=1,NSSO 00332 IF (JJSS-JPREV<1) THEN 00333 IL = IBOTTOM(JL) 00334 IPREV = JJSS-JPREV+NSSO 00335 ELSE 00336 IL = JL 00337 IPREV = JJSS-JPREV 00338 END IF 00339 ! no left point 00340 IF (IL==0) EXIT 00341 ! subgrid data found 00342 IF (LSSQO(JISS,IPREV,IL)) EXIT 00343 END DO 00344 00345 IF (.NOT. (JPREV>=NSSO+1 .OR. IL==0)) THEN 00346 ZDZSDY(JISS,1:JJSS) = ( XSSQO(JISS,JJSS,JL) - XSSQO(JISS,IPREV,IL)) & 00347 / FLOAT(JPREV) / ZDYEFF 00348 ! 00349 GDZSDY(JISS,1:JJSS) = .TRUE. 00350 END IF 00351 END IF 00352 ! 00353 ! 00354 GBOUND=.FALSE. 00355 ! 00356 ! 00357 ! 00358 !* 4.6 end of loop on jss index 00359 ! ------------------------ 00360 ! 00361 END DO 00362 ! 00363 !* 4.7 end of loop on iss index 00364 ! ------------------------ 00365 ! 00366 END DO 00367 !---------------------------------------------------------------------------- 00368 ! 00369 !* 5. Computations of tensor terms 00370 ! ---------------------------- 00371 ! 00372 ! 00373 !* 5.1 test to know if term Hxy is computable 00374 ! -------------------------------------- 00375 ! 00376 !* 2 values are necessary in the grid point to compute anisotropy 00377 ! 00378 IF ( COUNT(GDZSDX(:,:).AND.GDZSDY(:,:)) ==0 ) CYCLE 00379 ! 00380 ! 00381 !* 5.2 SSO quantities are computable 00382 ! ----------------------------- 00383 00384 OSSO(JL)=.TRUE. 00385 OSSO_ANIS(JL)=COUNT(GDZSDX(:,:).AND.GDZSDY(:,:))>1 00386 ! 00387 ! 00388 !* 5.3 term Hxx 00389 ! -------- 00390 ! 00391 ZHXX(JL) = SUM(ZDZSDX(:,:)*ZDZSDX(:,:),MASK=GDZSDX(:,:).AND.GDZSDY(:,:))& 00392 /COUNT(GDZSDX(:,:).AND.GDZSDY(:,:)) 00393 ! 00394 !* 5.4 term Hyy 00395 ! -------- 00396 ! 00397 ZHYY(JL) = SUM(ZDZSDY(:,:)*ZDZSDY(:,:),MASK=GDZSDX(:,:).AND.GDZSDY(:,:))& 00398 /COUNT(GDZSDX(:,:).AND.GDZSDY(:,:)) 00399 ! 00400 !* 5.5 term Hxy 00401 ! -------- 00402 ! 00403 ZHXY(JL) = SUM(ZDZSDX(:,:)*ZDZSDY(:,:),MASK=GDZSDX(:,:).AND.GDZSDY(:,:))& 00404 /COUNT(GDZSDX(:,:).AND.GDZSDY(:,:)) 00405 ! 00406 !------------------------------------------------------------------------------- 00407 ! 00408 !* 6. Next MESONH grid point 00409 ! ---------------------- 00410 ! 00411 END DO 00412 ! 00413 !------------------------------------------------------------------------------- 00414 ! 00415 !* 7. Diagonalization of the tensor 00416 ! ----------------------------- 00417 ! 00418 WHERE (OSSO(:)) 00419 ZK(:)=0.5*(ZHXX(:)+ZHYY(:)) 00420 ZL(:)=0.5*(ZHXX(:)-ZHYY(:)) 00421 ZM(:)= ZHXY(:) 00422 END WHERE 00423 ! 00424 !------------------------------------------------------------------------------- 00425 ! 00426 !* 8. S.S.O. characteristics 00427 ! ---------------------- 00428 ! 00429 !* 8.1 S.S.O. direction of main axis 00430 ! ----------------------------- 00431 ! 00432 WHERE (OSSO(:) .AND. ZL(:)>1.E-30 ) 00433 XSSO_DIR(:) = 0.5* ATAN(ZM/ZL) * (180./XPI) 00434 END WHERE 00435 ! 00436 WHERE (OSSO(:) .AND. ZL(:)<-1.E-30 ) 00437 XSSO_DIR(:) = 0.5* ATAN(ZM/ZL) * (180./XPI) + 90. 00438 END WHERE 00439 ! 00440 WHERE (OSSO(:) .AND. ABS(ZL(:))<=1.E-30 ) 00441 XSSO_DIR(:) = 45. 00442 END WHERE 00443 ! 00444 WHERE (OSSO(:) .AND. XSSO_DIR(:)>90. ) 00445 XSSO_DIR(:) = XSSO_DIR(:) - 180. 00446 END WHERE 00447 ! 00448 !* 8.2 S.S.O. slope 00449 ! ------------ 00450 ! 00451 WHERE (OSSO(:)) 00452 XSSO_SLOPE(:) = SQRT( ZK+SQRT(ZL*ZL+ZM*ZM) ) 00453 END WHERE 00454 ! 00455 !* 8.3 S.S.O. anisotropy 00456 ! ----------------- 00457 ! 00458 WHERE (OSSO_ANIS(:) .AND. (ZK+SQRT(ZL*ZL+ZM*ZM)) >0. ) 00459 XSSO_ANIS(:)=SQRT( MAX(ZK-SQRT(ZL*ZL+ZM*ZM),0.) / (ZK+SQRT(ZL*ZL+ZM*ZM))) 00460 END WHERE 00461 ! 00462 WHERE (OSSO_ANIS(:) .AND. (ZK+SQRT(ZL*ZL+ZM*ZM))==0. ) 00463 XSSO_ANIS(:)=1. 00464 END WHERE 00465 IF (LHOOK) CALL DR_HOOK('SSO',1,ZHOOK_HANDLE) 00466 ! 00467 !------------------------------------------------------------------------------- 00468 ! 00469 END SUBROUTINE SSO
1.8.0
