New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
usrdef_zgr.F90 in NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC/usrdef_zgr.F90 @ 13562

Last change on this file since 13562 was 13562, checked in by gm, 4 years ago

zgr_pse created only for NS coast
File size: 26.1 KB
Line 
1MODULE usrdef_zgr
2   !!======================================================================
3   !!                       ***  MODULE  usrdef_zgr  ***
4   !!
5   !!                       ===  AM98 configuration  ===
6   !!
7   !! User defined : vertical coordinate system of a user configuration
8   !!======================================================================
9   !! History :  4.0  ! 2016-06  (G. Madec)  Original code
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   usr_def_zgr   : user defined vertical coordinate system
14   !!      zgr_z      : reference 1D z-coordinate
15   !!      zgr_top_bot: ocean top and bottom level indices
16   !!      zgr_zco    : 3D verticl coordinate in pure z-coordinate case
17   !!---------------------------------------------------------------------
18   USE oce            ! ocean variables
19   USE dom_oce        ! ocean domain
20   USE depth_e3       ! depth <=> e3
21   USE usrdef_nam
22   !
23   USE in_out_manager ! I/O manager
24   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
25   USE lib_mpp        ! distributed memory computing library
26
27   USE usrdef_nam , ONLY : rn_dx, nn_AM98, rn_abp, r1_abp    ! use penalisation parameter
28
29   IMPLICIT NONE
30   PRIVATE
31   !
32   INTEGER, PARAMETER, PRIVATE  ::   nT  = 1
33   INTEGER, PARAMETER, PRIVATE  ::   nU  = 2
34   INTEGER, PARAMETER, PRIVATE  ::   nV  = 3
35   INTEGER, PARAMETER, PRIVATE  ::   nF  = 4
36   !
37   PUBLIC   usr_def_zgr        ! called by domzgr.F90
38
39   !!----------------------------------------------------------------------
40   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
41   !! $Id: usrdef_zgr.F90 10425 2018-12-19 21:54:16Z smasson $
42   !! Software governed by the CeCILL license (see ./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
47      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
48      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
49      &                    pe3t  , pe3u  , pe3v   , pe3f ,             &   ! vertical scale factors
50      &                    pe3w  , pe3uw , pe3vw         ,             &   !     -      -      -
51      &                    k_top , k_bot                 )                 ! top & bottom ocean level
52      !!---------------------------------------------------------------------
53      !!              ***  ROUTINE usr_def_zgr  ***
54      !!
55      !! ** Purpose :   User defined the vertical coordinates
56      !!
57      !!----------------------------------------------------------------------
58      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
59      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
60      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
61      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
62      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
63      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
64      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors
65      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
66      !
67      INTEGER  ::   inum   ! local logical unit
68      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
69      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
70      !!----------------------------------------------------------------------
71      !
72      IF(lwp) WRITE(numout,*)
73      IF(lwp) WRITE(numout,*) 'usr_def_zgr : AM98 configuration (z-coordinate closed flat box ocean without cavities)'
74      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
75      !
76      !
77      ! type of vertical coordinate
78      ! ---------------------------
79      ld_zco    = .FALSE.         ! AM98 case:  z-coordinate without ocean cavities
80      ld_zps    = .FALSE.
81      ld_sco    = .TRUE.
82      ld_isfcav = .FALSE.
83      !
84      !
85      ! Build the vertical coordinate system
86      ! ------------------------------------
87      CALL zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
88      !
89      CALL zgr_msk_top_bot( k_top , k_bot)                 ! masked top and bottom ocean t-level indices
90      !
91      !                                                     ! z-coordinate (3D arrays) from the 1D z-coord.
92      CALL zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in  : 1D reference vertical coordinate
93         &          pdept   , pdepw   ,                     &   ! out : 3D t & w-points depth
94         &          pe3t    , pe3u    , pe3v   , pe3f   ,   &   !       vertical scale factors
95         &          pe3w    , pe3uw   , pe3vw             )     !           -      -      -
96      !
97   END SUBROUTINE usr_def_zgr
98
99
100   SUBROUTINE zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
101      !!----------------------------------------------------------------------
102      !!                   ***  ROUTINE zgr_z  ***
103      !!
104      !! ** Purpose :   set the 1D depth of model levels and the resulting
105      !!              vertical scale factors.
106      !!
107      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
108      !!       The depth of model levels is set from dep(k), an analytical function:
109      !!                   w-level: depw_1d  = dep(k)
110      !!                   t-level: dept_1d  = dep(k+0.5)
111      !!       The scale factors are the discrete derivative of the depth:
112      !!                   e3w_1d(jk) = dk[ dept_1d ]
113      !!                   e3t_1d(jk) = dk[ depw_1d ]
114      !!           with at top and bottom :
115      !!                   e3w_1d( 1 ) = 2 * ( dept_1d( 1 ) - depw_1d( 1 ) )
116      !!                   e3t_1d(jpk) = 2 * ( dept_1d(jpk) - depw_1d(jpk) )
117      !!       The depth are then re-computed from the sum of e3. This ensures
118      !!    that depths are identical when reading domain configuration file.
119      !!    Indeed, only e3. are saved in this file, depth are compute by a call
120      !!    to the e3_to_depth subroutine.
121      !!
122      !!       Here the Madec & Imbard (1996) function is used.
123      !!
124      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
125      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
126      !!
127      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
128      !!             Madec and Imbard, 1996, Clim. Dyn.
129      !!----------------------------------------------------------------------
130      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
131      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
132      !
133      INTEGER  ::   jk       ! dummy loop indices
134        !!----------------------------------------------------------------------
135      !
136        !
137      IF(lwp) THEN            ! Parameter print
138         WRITE(numout,*)
139         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates '
140         WRITE(numout,*) '    ~~~~~~~'
141      ENDIF
142
143      !
144      ! 1D Reference z-coordinate    (using Madec & Imbard 1996 function)
145      ! -------------------------
146      !
147      ! depth at T and W-points
148      pdepw_1d(1) =   0._wp
149      pdept_1d(1) = 250._wp
150      !
151      pdepw_1d(2) = 500._wp
152      pdept_1d(2) = 750._wp
153      !
154      !                       ! e3t and e3w from depth
155      CALL depth_to_e3( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d )
156      !
157      !                       ! recompute depths from SUM(e3)  <== needed
158      CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )
159      !
160      IF(lwp) THEN                        ! control print
161         WRITE(numout,*)
162         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
163         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
164         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
165      ENDIF
166      !
167   END SUBROUTINE zgr_z
168
169
170
171
172   SUBROUTINE zgr_msk_top_bot( k_top , k_bot)
173      !!----------------------------------------------------------------------
174      !!                    ***  ROUTINE zgr_msk_top_bot  ***
175      !!
176      !! ** Purpose :   set the masked top and bottom ocean t-levels
177      !!
178      !! ** Method  :   AM98 case = closed flat box ocean without ocean cavities
179      !!                   k_top = 1     except along north, south, east and west boundaries
180      !!                   k_bot = jpk-1 except along north, south, east and west boundaries
181      !!
182      !! ** Action  : - k_top : first wet ocean level index
183      !!              - k_bot : last  wet ocean level index
184      !!----------------------------------------------------------------------
185      INTEGER , DIMENSION(:,:), INTENT(out) ::   k_top , k_bot   ! first & last wet ocean level
186      !
187      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D local workspace
188      INTEGER  ::   ji, jj,jk, jl                    ! dummy loop indices
189      REAL(wp) ::   zylim0, zylim1, zxlim0, zxlim1 ! limit of the domain [m]
190      REAL(WP) ::   zcoeff    ! local scalar
191      !!----------------------------------------------------------------------
192      !
193      IF(lwp) WRITE(numout,*)
194      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom wet ocean levels.'
195      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
196      IF(lwp) WRITE(numout,*) '       AM98 case : closed flat box ocean without ocean cavities'
197      !
198      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
199      !
200      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed)
201      !
202      !
203# if defined key_bvp
204      ! Faire la pénalisation
205      ! ktop = where phi flotte < 1% puis envoie au mask
206      !  Penalize the first inner layer of fluid
207      ! -----------------------------------------------------
208      zylim0 =  10000._wp     ! +   10km
209      zylim1 = 2010000._wp    !   2010km
210      zxlim0 = -200000._wp    ! -  200km   - zgr_pse
211      zxlim1 = 2010000._wp    !   2010km
212      rpo (:,:,:) = rn_abp  ;   rpou(:,:,:) = rn_abp  ;   rpov(:,:,:) = rn_abp   ;   rpof(:,:,:) = rn_abp
213      !
214      DO jj = 2, jpjm1
215         DO ji = 2, jpim1
216           ! it is "largely" whithin the box
217           IF ( gphit(ji,jj) > zylim0 .AND. gphit(ji,jj) < zylim1 .AND. &
218              & glamt(ji,jj) > zxlim0 .AND. glamt(ji,jj) < zxlim1       )  THEN
219              CALL zgr_pse (ji,jj,1,glamf,gphif,rpo , nT )
220              CALL zgr_pse (ji,jj,1,glamt,gphit,rpof, nF )
221              CALL zgr_pse (ji,jj,1,glamv,gphiv,rpou, nU )
222              CALL zgr_pse (ji,jj,1,glamu,gphiu,rpov, nV )
223           END IF
224        END DO
225      END DO
226      !
227      r1_rpo (:,:,:) = r1_abp ; r1_rpou(:,:,:) = r1_abp ; r1_rpov(:,:,:) = r1_abp
228      bmpt   (:,:,:) = 0._wp  ; bmpu   (:,:,:) = 0._wp  ; bmpv   (:,:,:) = 0._wp
229      !
230      r1_rpof(:,:,:) = r1_abp
231      !
232      DO jj = 1, jpj
233         DO ji = 1, jpi
234           !
235           ! thinner
236           !bmpu(ji,jj, 1:jpkm1) = rn_fsp * ( 1._wp - rpou(ji,jj,1:jpkm1) )
237           !bmpv(ji,jj, 1:jpkm1) = rn_fsp * ( 1._wp - rpov(ji,jj,1:jpkm1) )
238           ! classic constant
239           IF( rpou(ji,jj,1) < 1._wp )   bmpu(ji,jj, 1:jpkm1) = rn_fsp
240           IF( rpov(ji,jj,1) < 1._wp )   bmpv(ji,jj, 1:jpkm1) = rn_fsp
241           !
242           r1_rpo (ji,jj,1:jpkm1) = 1._wp / rpo (ji,jj,1:jpkm1)
243           r1_rpou(ji,jj,1:jpkm1) = 1._wp / rpou(ji,jj,1:jpkm1)
244           r1_rpov(ji,jj,1:jpkm1) = 1._wp / rpov(ji,jj,1:jpkm1)
245           r1_rpof(ji,jj,1:jpkm1) = 1._wp / rpof(ji,jj,1:jpkm1)
246        END DO
247      END DO
248      !
249      WHERE(rpo(:,:,1) <= rn_abp )
250        rpo  (:,:,1) = rn_abp
251      END WHERE
252      WHERE(rpou(:,:,1) <= rn_abp )
253        rpou  (:,:,1) = rn_abp
254      END WHERE
255      WHERE(rpov(:,:,1) <= rn_abp )
256        rpov  (:,:,1) = rn_abp
257      END WHERE
258      WHERE(rpof(:,:,1) <= rn_abp )
259        rpof  (:,:,1) = rn_abp
260      END WHERE
261      !
262      ! definition du mask
263      k_top(:,:) = 1    ! = ocean
264      k_bot(:,:) = NINT( z2d(:,:) )
265      WHERE(rpo(:,:,1) <= rn_abp )
266        k_top(:,:) = 0    ! = land
267        k_bot(:,:) = 0
268      END WHERE
269      !
270      CALL lbc_lnk_multi( 'usrdef_zgr', rpo , 'T', 1._wp, r1_rpo  , 'T', 1._wp,                   &
271             &                          rpou, 'U', 1._wp, r1_rpou , 'U', 1._wp, bmpu, 'U', 1._wp,    &
272             &                          rpov, 'V', 1._wp, r1_rpov , 'V', 1._wp, bmpv, 'V', 1._wp,    &
273             &                          rpof, 'F', 1._wp, r1_rpof , 'F', 1._wp, kfillmode=jpfillcopy )
274# else
275      ! Dans l'idéal, j'aimerais bien pouvoir me passer de cette boucle
276      zylim0 =   10000._wp    !  +10km
277      zylim1 = 2010000._wp    ! 2010km
278      zxlim0 =   10000._wp    !  +10km
279      zxlim1 = 2010000._wp    ! 2010km
280      !
281      DO jj = 1, jpj
282         DO ji = 1, jpi
283         ! if T point in the 2000 km x 2000 km domain
284         ! IF ( gphit(ji,jj) > zylim0 .AND. gphit(ji,jj) < zylim1 .AND. &
285         !   & glamt(ji,jj) > zxlim0 .AND. glamt(ji,jj) < zxlim1       )  THEN
286         ! if U,V points are in the 2000 km x 2000 km domain
287         IF ( gphiv(ji,jj) > zylim0 .AND. gphiv(ji,jj) < zylim1 .AND. &
288            & glamu(ji,jj) > zxlim0 .AND. glamu(ji,jj) < zxlim1       )  THEN
289         k_top(ji,jj) = 1    ! = ocean
290         k_bot(ji,jj) = NINT( z2d(ji,jj) )
291         ELSE
292         k_top(ji,jj) = 0    ! = land
293         k_bot(ji,jj) = 0
294         END IF
295         END DO
296      END DO
297      ! mask the lonely corners
298      DO jj = 2, jpim1
299         DO ji = 2, jpim1
300         zcoeff = k_top(ji+1,jj) + k_top(ji,jj+1)   &
301            +     k_top(ji-1,jj) + k_top(ji,jj-1)
302         IF ( zcoeff <= 1._wp )   THEN
303            k_top(ji,jj) = 0    ! = land
304            k_bot(ji,jj) = 0
305         END IF
306         END DO
307      END DO
308
309# endif
310
311   END SUBROUTINE zgr_msk_top_bot
312
313
314   SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate
315      &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth
316      &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors
317      &                pe3w    , pe3uw   , pe3vw             )     !          -      -      -
318      !!----------------------------------------------------------------------
319      !!                  ***  ROUTINE zgr_zco  ***
320      !!
321      !! ** Purpose :   define the reference z-coordinate system
322      !!
323      !! ** Method  :   set 3D coord. arrays to reference 1D array
324      !!----------------------------------------------------------------------
325      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
326      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
327      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m]
328      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
329      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
330      !
331      INTEGER  ::   jk
332      !!----------------------------------------------------------------------
333      !
334      DO jk = 1, jpk
335         pdept(:,:,jk) = pdept_1d(jk)
336         pdepw(:,:,jk) = pdepw_1d(jk)
337         pe3t (:,:,jk) = pe3t_1d (jk)
338         pe3u (:,:,jk) = pe3t_1d (jk)
339         pe3v (:,:,jk) = pe3t_1d (jk)
340         pe3f (:,:,jk) = pe3t_1d (jk)
341         pe3w (:,:,jk) = pe3w_1d (jk)
342         pe3uw(:,:,jk) = pe3w_1d (jk)
343         pe3vw(:,:,jk) = pe3w_1d (jk)
344      END DO
345      !
346   END SUBROUTINE zgr_zco
347
348
349   SUBROUTINE zgr_pse( ki, kj, kk,                  &
350      &                pglam,pgphi,prpo,            &
351      &                cpoint                       )
352      !!----------------------------------------------------------------------
353      !!                  ***  ROUTINE zgr_pse  ***
354      !!
355      !! ** Purpose :   Estimate the porosity field within the cell
356      !!
357      !! ** Method  :   find the intersection points with the boundaries
358      !!                calculate the area land/water
359      !!----------------------------------------------------------------------
360      INTEGER                   , INTENT(in   ) ::   ki, kj, kk    ! coordinate of the center of the cell
361      REAL, DIMENSION(:,:)      , INTENT(in   ) ::   pglam, pgphi  ! grid-point position          [m]
362      REAL, DIMENSION(:,:,:)    , INTENT(inout) ::   prpo          ! porosity field
363      INTEGER                   , INTENT(in   ) ::   cpoint        ! type of point in the C-grid
364                                                                   ! T : 1, U : 2, V : 3, F : 4
365      INTEGER  ::  jtype, jlen                     ! geometrical figure ; nearest land(=0)/water(=1) point
366      REAL ::   z1d, z1_dx, zl                     ! dummy variable
367      REAL, DIMENSION(2) ::   zA, zB, zC, zD, zM   ! coordinate array M(x,y)
368      !
369      INTEGER  ::   jkM                      ! local index of intersection points
370      REAL,    DIMENSION(2)    :: zlm        ! save the lenght alond the edges
371      INTEGER, DIMENSION(4)    :: zsp, zss   ! Sign Points, sum sign
372                                             ! define the positions of points to the limit
373      !!----------------------------------------------------------------------
374      !
375      z1_dx =   REAL(nn_AM98, wp) / rn_dx       ! [1/m] inverse gridspacing used
376      !REAL :: dx = 1.
377      !REAL :: z1_dx =   1.       ! [1/m] inverse gridspacing used
378      !
379      ! WRITE(*,*) 'BEGINNING of the routine zgr_pse'
380      !                                      !==  Preparatory work  ==!
381      SELECT CASE ( cpoint )                     !* Defining vertices
382        CASE ( nT )                                               ! F coordinates
383          zA(1) = pglam(ki-1,kj-1 ) ; zA(2) = pgphi(ki-1,kj-1 )
384          zB(1) = pglam(ki  ,kj-1 ) ; zB(2) = pgphi(ki  ,kj-1 )
385          zC(1) = pglam(ki  ,kj   ) ; zC(2) = pgphi(ki  ,kj   )
386          zD(1) = pglam(ki-1,kj   ) ; zD(2) = pgphi(ki-1,kj   )
387        CASE ( nF )                                               ! T coordinates
388          zA(1) = pglam(ki  ,kj   ) ; zA(2) = pgphi(ki  ,kj   )
389          zB(1) = pglam(ki+1,kj   ) ; zB(2) = pgphi(ki+1,kj   )
390          zC(1) = pglam(ki+1,kj+1 ) ; zC(2) = pgphi(ki+1,kj+1 )
391          zD(1) = pglam(ki  ,kj+1 ) ; zD(2) = pgphi(ki  ,kj+1 )
392        CASE ( nU )                                               ! V coordinates
393          zA(1) = pglam(ki  ,kj-1 ) ; zA(2) = pgphi(ki  ,kj-1 )
394          zB(1) = pglam(ki+1,kj-1 ) ; zB(2) = pgphi(ki+1,kj-1 )
395          zC(1) = pglam(ki+1,kj   ) ; zC(2) = pgphi(ki+1,kj   )
396          zD(1) = pglam(ki  ,kj   ) ; zD(2) = pgphi(ki  ,kj   )
397        CASE ( nV )                                               ! U coordinates
398          zA(1) = pglam(ki-1,kj   ) ; zA(2) = pgphi(ki-1,kj   )
399          zB(1) = pglam(ki  ,kj   ) ; zB(2) = pgphi(ki  ,kj   )
400          zC(1) = pglam(ki  ,kj+1 ) ; zC(2) = pgphi(ki  ,kj+1 )
401          zD(1) = pglam(ki-1,kj+1 ) ; zD(2) = pgphi(ki-1,kj+1 )
402      END SELECT
403      !
404        ! WRITE(*,*) ' zA = ',zA
405        ! WRITE(*,*) ' zB = ',zB
406        ! WRITE(*,*) ' zC = ',zC
407        ! WRITE(*,*) ' zD = ',zD
408      !
409      !                                       zss(3)
410      !                                zsp(4)   |     zsp(3)
411      !      D --------- C                 D --------- C
412      !      |           |                 |           |
413      !      |     +     |         zss(4)--|     +     |-- zss(2)
414      !      |  (ki,kj)  |                 |           |
415      !      A --------- B                 A --------- B
416      !                                zsp(1)   |     zsp(2)
417      !                                       zss(1)
418      !
419      ! !! case 1 : (x=)glam = 0             !* Defining position to the boundary
420      ! !! 1 = is in water ; 0 in land
421      ! zsp(:) = 0   ;   zss(:) = 0
422      ! IF ( zA(1) > ( 0. )  )   zsp(1) = 1    ! A
423      ! IF ( zB(1) > ( 0. )  )   zsp(2) = 1    ! B
424      ! IF ( zC(1) > ( 0. )  )   zsp(3) = 1    ! C
425      ! IF ( zD(1) > ( 0. )  )   zsp(4) = 1    ! D
426      ! !
427      !! case 1 : (x=)glam = 0             !* Defining position to the boundary
428      !! 1 = is in water ; 0 in land
429      zsp(:) = 0   ;   zss(:) = 0
430      IF ( zA(1) > ( 0. )  )   zsp(1) = 1    ! A
431      IF ( zB(1) > ( 0. )  )   zsp(2) = 1    ! B
432      IF ( zC(1) > ( 0. )  )   zsp(3) = 1    ! C
433      IF ( zD(1) > ( 0. )  )   zsp(4) = 1    ! D
434      !
435      ! !! case 2 : (x=)gphi = 2000 km       !* Defining position to the boundary
436      ! !! 1 = is in water ; 0 in land
437      ! zsp(:) = 0   ;   zss(:) = 0
438      ! IF ( zA(2) < ( 200000. )  )   zsp(1) = 1    ! A
439      ! IF ( zB(2) < ( 200000. )  )   zsp(2) = 1    ! B
440      ! IF ( zC(2) < ( 200000. )  )   zsp(3) = 1    ! C
441      ! IF ( zD(2) < ( 200000. )  )   zsp(4) = 1    ! D
442      !
443      zss(1) = zsp(1) + zsp(2)
444      zss(2) = zsp(2) + zsp(3)
445      zss(3) = zsp(3) + zsp(4)
446      zss(4) = zsp(4) + zsp(1)
447      !
448      IF      ( SUM( zsp(:) ) == 3 )  THEN
449         jtype = 1   ;   jlen = 0
450         ! WRITE(*,*) 'land triangle'
451      ELSE IF ( SUM( zsp(:) ) == 2 )  THEN
452         jtype = 2   ;  jlen = 1
453         ! WRITE(*,*) 'water trapeze'
454      ELSE IF ( SUM( zsp(:) ) == 1 )  THEN
455         jtype = 3    ;  jlen = 1
456         ! WRITE(*,*) 'water triangle'
457      ELSE IF ( pglam(ki,kj) > 0._wp ) THEN
458         ! SUM( zsp(:) ) == 4
459         ! prpo(ki,kj,kk) = prpo(ki,kj,kk) + 1._wp
460         prpo(ki,kj,kk) =  1._wp
461         ! full water
462         RETURN
463      ELSE
464        ! full land
465        RETURN
466      ENDIF
467      !
468      ! WRITE(*,*) ' zsp = ',zsp(4), zsp(3)
469      ! WRITE(*,*) '       ',zsp(1), zsp(2)
470      !
471      zlm(:) = 0.
472      !                            !==  Searching for the intersection points  ==!
473      !!                                !* NS west boundary
474      zM(1) = 0.
475      ! WRITE(*,*) ' zlm = ',zlm
476      jkM = 1
477       !    A ------- B
478      IF ( zss(1) == 1 ) THEN
479        zM(2) = weightedMx(zA,zB, z1_dx, 0.)
480        !
481        IF (zsp(1) == jlen) THEN
482          ! WRITE(*,*) ' AM '
483          zlm(jkM) = lenght(zA,zM)
484        ELSE
485          ! WRITE(*,*) ' MB '
486          zlm(jkM) = lenght(zB,zM)
487        ENDIF
488        jkM = jkM + 1
489        ! WRITE(*,*) ' '
490        ! WRITE(*,*) ' zM = ', zM
491        ! WRITE(*,*) ' zlm = ',zlm
492      ENDIF
493       !   C
494       !   |
495       !   |
496       !   |
497       !   B
498      IF ( zss(2) == 1 ) THEN
499        zM(2) = weightedMx(zB,zC, z1_dx,0.)
500        !
501        IF (zsp(2) == jlen) THEN
502          zlm(jkM) = lenght(zB,zM)
503        ELSE
504          zlm(jkM) = lenght(zC,zM)
505        ENDIF
506        jkM = jkM + 1
507        ! WRITE(*,*) ' '
508        ! WRITE(*,*) ' BC '
509        ! WRITE(*,*) ' zM = ', zM
510        ! WRITE(*,*) ' zlm = ',zlm
511      ENDIF
512      !  !      D ------- C
513      IF ( zss(3) == 1 ) THEN
514        zM(2) = weightedMx(zC,zD, z1_dx,0.)
515        !
516        IF (zsp(3) == jlen) THEN
517          ! WRITE(*,*) ' CM '
518          zlm(jkM) = lenght(zC,zM)
519        ELSE
520          ! WRITE(*,*) ' MD '
521          zlm(jkM) = lenght(zD,zM)
522        ENDIF
523        jkM = jkM + 1
524        ! WRITE(*,*) ' '
525        ! WRITE(*,*) ' zM = ', zM
526        ! WRITE(*,*) ' zlm = ',zlm
527      ENDIF
528      !  !   D
529      !      |
530      !      |
531      !      |
532      !      A
533      IF ( zss(4) == 1 ) THEN
534        zM(2) = weightedMx(zD,zA, z1_dx, 0.)
535        !
536        IF (zsp(4) == jlen) THEN
537          zlm(jkM) = lenght(zD,zM)
538        ELSE
539          zlm(jkM) = lenght(zA,zM)
540        ENDIF
541        jkM = jkM + 1
542        ! WRITE(*,*) ' '
543        ! WRITE(*,*) ' DA '
544        ! WRITE(*,*) ' zM = ', zM
545        ! WRITE(*,*) ' zlm = ',zlm
546      ENDIF
547      ! WRITE(*,*) 'ENDIF  zlm = ',zlm
548      !
549      !
550      !                           !==  Calculating land area  ==!
551      !  ratio of water upon the surface
552      IF     (jtype == 1) THEN    ! land triangle
553        z1d = 1. - 0.5 * zlm(1) * zlm(2) * z1_dx * z1_dx
554      ELSEIF (jtype == 2) THEN    ! water trapeze
555        z1d =      0.5 * ( zlm(1) + zlm(2) ) * z1_dx
556      ELSEIF (jtype == 3) THEN    ! water triangle
557        z1d =      0.5 * zlm(1) * zlm(2) * z1_dx * z1_dx
558      ENDIF
559      ! update the porosity field in += manner
560      ! prpo(ki,kj,kk) = prpo(ki,kj,kk) + z1d
561      prpo(ki,kj,kk) = z1d
562      ! WRITE(*,*) 'prpo(ki,kj,kk)=',prpo(ki,kj,kk)
563      ! ...
564      ! NS east boundary
565      ! WE south boundary
566      ! WE north boundary
567      !
568      !
569      ! WRITE(*,*) 'END of the routine zgr_pse'
570      !
571   END SUBROUTINE zgr_pse
572
573
574   FUNCTION weightedMx(A,B,z1_dx,xl)  RESULT(f)
575      !!----------------------------------------------------------------------
576      !!                 ***  ROUTINE weightedM  ***
577      !!
578      !! ** Purpose :
579      !!
580      !! ** Method  :
581      !!
582      !!----------------------------------------------------------------------
583      IMPLICIT NONE
584      REAL, DIMENSION(2), INTENT(in) :: A,B
585      REAL,                 INTENT(in) :: z1_dx, xl
586      REAL                             ::  f                 ! result
587      !!----------------------------------------------------------------------
588      !
589      f  =    (   A(2) * ( B(1) - xl     )    &
590         &    +   B(2) * ( xl     - A(1) )    )   &
591         &    /            ( B(1) - A(1) )
592      !
593   END FUNCTION weightedMx
594
595   FUNCTION lenght(A,B)  RESULT(f)
596      !!----------------------------------------------------------------------
597      !!                 ***  ROUTINE weighted_lenght_X  ***
598      !!
599      !! ** Purpose :
600      !!
601      !! ** Method  :
602      !!
603      !!----------------------------------------------------------------------
604      IMPLICIT NONE
605      REAL, DIMENSION(2), INTENT(in) :: A,B
606      REAL             ::  f                   ! result
607      !!----------------------------------------------------------------------
608      !
609      f  = SQRT(   ( A(1) - B(1) ) * ( A(1) - B(1) )    &
610         &     +   ( A(2) - B(2) ) * ( A(2) - B(2) )    )
611      !
612   END FUNCTION lenght
613
614
615   !!======================================================================
616END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.