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.
Changeset 7698 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 – NEMO

Ignore:
Timestamp:
2017-02-18T10:02:03+01:00 (7 years ago)
Author:
mocavero
Message:

update trunk with OpenMP parallelization

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r7646 r7698  
    115115      REAL(wp), POINTER, DIMENSION(:,:)   ::   opning          ! rate of opening due to divergence/shear 
    116116      REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
     117      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_ai 
    117118      ! 
    118119      INTEGER, PARAMETER ::   nitermax = 20     
    … …  
    122123      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    123124 
    124       CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross ) 
     125      CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, z_ai ) 
    125126 
    126127      ! conservation test 
    … …  
    135136      ! 
    136137 
     138!$OMP PARALLEL DO schedule(static) private(jj,ji) 
    137139      DO jj = 1, jpj                                     ! Initialize arrays. 
    138140         DO ji = 1, jpi 
    … …  
    192194         !  closing rate to a gross closing rate.   
    193195         ! NOTE: 0 < aksum <= 1 
    194          closing_gross(:,:) = closing_net(:,:) / aksum(:,:) 
     196!$OMP PARALLEL 
     197!$OMP DO schedule(static) private(jj,ji) 
     198         DO jj = 1, jpj 
     199            DO ji = 1, jpi 
     200               closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 
     201            END DO 
     202         END DO 
    195203 
    196204         ! correction to closing rate and opening if closing rate is excessive 
    … …  
    198206         ! Reduce the closing rate if more than 100% of the open water  
    199207         ! would be removed.  Reduce the opening rate proportionately. 
     208!$OMP DO schedule(static) private(jj,ji,za,zfac) 
    200209         DO jj = 1, jpj 
    201210            DO ji = 1, jpi 
    … …  
    216225         ! would be removed.  Reduce the opening rate proportionately. 
    217226         DO jl = 1, jpl 
     227!$OMP DO schedule(static) private(jj,ji,za,zfac) 
    218228            DO jj = 1, jpj 
    219229               DO ji = 1, jpi 
    … …  
    226236            END DO 
    227237         END DO 
     238!$OMP END PARALLEL 
    228239 
    229240         ! 3.3 Redistribute area, volume, and energy. 
    … …  
    236247         !-----------------------------------------------------------------------------! 
    237248         ! This is in general not equal to one because of divergence during transport 
    238          asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
     249!$OMP PARALLEL 
     250!$OMP DO schedule(static) private(jj,ji) 
     251         DO jj = 1, jpj 
     252            DO ji = 1, jpi 
     253               asum(ji,jj) = 0._wp 
     254               z_ai(ji,jj) = 0._wp 
     255            END DO 
     256         END DO 
     257         DO jl = 1, jpl 
     258!$OMP DO schedule(static) private(jj,ji) 
     259            DO jj = 1, jpj 
     260               DO ji = 1, jpi 
     261                  z_ai(ji,jj) = z_ai(ji,jj) + a_i(ji,jj,jl) 
     262               END DO 
     263            END DO 
     264         END DO 
     265!$OMP DO schedule(static) private(jj,ji) 
     266         DO jj = 1, jpj 
     267            DO ji = 1, jpi 
     268               asum(ji,jj) = ato_i(ji,jj) + z_ai(ji,jj) 
     269            END DO 
     270         END DO 
    239271 
    240272         ! 3.5 Do we keep on iterating ??? 
    … …  
    244276 
    245277         iterate_ridging = 0 
     278!$OMP DO schedule(static) private(jj,ji) 
    246279         DO jj = 1, jpj 
    247280            DO ji = 1, jpi 
    … …  
    258291            END DO 
    259292         END DO 
     293!$OMP END PARALLEL 
    260294 
    261295         IF( lk_mpp )   CALL mpp_max( iterate_ridging ) 
    … …  
    289323      IF( ln_ctl )       CALL lim_prt3D( 'limitd_me' ) 
    290324 
    291       CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross ) 
     325      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, z_ai ) 
    292326      ! 
    293327      IF( nn_timing == 1 )  CALL timing_stop('limitd_me') 
    … …  
    306340      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar 
    307341      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n 
     342      REAL(wp), POINTER, DIMENSION(:,:) ::   z_ai 
    308343      !------------------------------------------------------------------------------! 
    309344 
    310345      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 
     346      CALL wrk_alloc( jpi,jpj,z_ai ) 
    311347 
    312348      Gstari     = 1.0/rn_gstar     
    313349      astari     = 1.0/rn_astar     
    314       aksum(:,:)    = 0.0 
    315       athorn(:,:,:) = 0.0 
    316       aridge(:,:,:) = 0.0 
    317       araft (:,:,:) = 0.0 
     350!$OMP PARALLEL 
     351!$OMP DO schedule(static) private(jj,ji) 
     352      DO jj = 1, jpj 
     353         DO ji = 1, jpi 
     354            aksum(ji,jj)    = 0.0 
     355         END DO 
     356      END DO 
     357!$OMP END DO NOWAIT 
     358      DO jl = 1, jpl 
     359!$OMP DO schedule(static) private(jj,ji) 
     360         DO jj = 1, jpj 
     361            DO ji = 1, jpi 
     362               athorn(ji,jj,jl) = 0.0 
     363               aridge(ji,jj,jl) = 0.0 
     364               araft (ji,jj,jl) = 0.0 
     365            END DO 
     366         END DO 
     367      END DO 
     368!$OMP END PARALLEL 
    318369 
    319370      ! Zero out categories with very small areas 
    320371      CALL lim_var_zapsmall 
    321372 
     373!$OMP PARALLEL 
    322374      ! Ice thickness needed for rafting 
    323375      DO jl = 1, jpl 
     376!$OMP DO schedule(static) private(jj,ji,rswitch) 
    324377         DO jj = 1, jpj 
    325378            DO ji = 1, jpi 
    … …  
    336389      ! Compute total area of ice plus open water. 
    337390      ! This is in general not equal to one because of divergence during transport 
    338       asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
    339  
     391 
     392!$OMP DO schedule(static) private(jj,ji) 
     393         DO jj = 1, jpj 
     394            DO ji = 1, jpi 
     395               asum(ji,jj) = 0._wp 
     396               z_ai(ji,jj) = 0._wp 
     397            END DO 
     398         END DO 
     399         DO jl = 1, jpl 
     400!$OMP DO schedule(static) private(jj,ji) 
     401            DO jj = 1, jpj 
     402               DO ji = 1, jpi 
     403                  z_ai(ji,jj) = z_ai(ji,jj) + a_i(ji,jj,jl) 
     404               END DO 
     405            END DO 
     406         END DO 
     407!$OMP DO schedule(static) private(jj,ji) 
     408         DO jj = 1, jpj 
     409            DO ji = 1, jpi 
     410               asum(ji,jj) = ato_i(ji,jj) + z_ai(ji,jj) 
     411            END DO 
     412         END DO 
    340413      ! Compute cumulative thickness distribution function 
    341414      ! Compute the cumulative thickness distribution function Gsum, 
    342415      ! where Gsum(n) is the fractional area in categories 0 to n. 
    343416      ! initial value (in h = 0) equals open water area 
    344       Gsum(:,:,-1) = 0._wp 
    345       Gsum(:,:,0 ) = ato_i(:,:) 
     417!$OMP DO schedule(static) private(jj,ji) 
     418      DO jj = 1, jpj 
     419         DO ji = 1, jpi 
     420            Gsum(ji,jj,-1) = 0._wp 
     421            Gsum(ji,jj,0 ) = ato_i(ji,jj) 
     422         END DO 
     423      END DO 
    346424      ! for each value of h, you have to add ice concentration then 
    347425      DO jl = 1, jpl 
    348          Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 
     426!$OMP DO schedule(static) private(jj,ji) 
     427         DO jj = 1, jpj 
     428            DO ji = 1, jpi 
     429               Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
     430            END DO 
     431         END DO 
    349432      END DO 
    350433 
    351434      ! Normalize the cumulative distribution to 1 
    352435      DO jl = 0, jpl 
    353          Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:) 
     436!$OMP DO schedule(static) private(jj,ji) 
     437         DO jj = 1, jpj 
     438            DO ji = 1, jpi 
     439               Gsum(ji,jj,jl) = Gsum(ji,jj,jl) / asum(ji,jj) 
     440            END DO 
     441         END DO 
    354442      END DO 
     443!$OMP END PARALLEL 
    355444 
    356445      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) 
    … …  
    369458      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
    370459         DO jl = 0, jpl     
     460!$OMP PARALLEL DO schedule(static) private(jj,ji) 
    371461            DO jj = 1, jpj  
    372462               DO ji = 1, jpi 
    … …  
    387477         !                         
    388478         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array 
     479!$OMP PARALLEL 
    389480         DO jl = -1, jpl 
    390             Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 
     481!$OMP DO schedule(static) private(jj,ji) 
     482            DO jj = 1, jpj  
     483               DO ji = 1, jpi 
     484                  Gsum(ji,jj,jl) = EXP( -Gsum(ji,jj,jl) * astari ) * zdummy 
     485               END DO 
     486            END DO 
    391487         END DO 
    392488         DO jl = 0, jpl 
    393              athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 
    394          END DO 
     489!$OMP DO schedule(static) private(jj,ji) 
     490            DO jj = 1, jpj  
     491               DO ji = 1, jpi 
     492                  athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl) 
     493               END DO 
     494            END DO 
     495         END DO 
     496!$OMP END PARALLEL 
    395497         ! 
    396498      ENDIF 
    … …  
    400502         ! 
    401503         DO jl = 1, jpl 
     504!$OMP PARALLEL DO schedule(static) private(jj,ji,zdummy) 
    402505            DO jj = 1, jpj  
    403506               DO ji = 1, jpi 
    … …  
    412515         ! 
    413516         DO jl = 1, jpl 
    414             aridge(:,:,jl) = athorn(:,:,jl) 
     517!$OMP PARALLEL DO schedule(static) private(jj,ji) 
     518            DO jj = 1, jpj  
     519               DO ji = 1, jpi 
     520                  aridge(ji,jj,jl) = athorn(ji,jj,jl) 
     521               END DO 
     522            END DO 
    415523         END DO 
    416524         ! 
    … …  
    418526         ! 
    419527         DO jl = 1, jpl 
    420             araft(:,:,jl) = athorn(:,:,jl) 
     528!$OMP PARALLEL DO schedule(static) private(jj,ji) 
     529            DO jj = 1, jpj  
     530               DO ji = 1, jpi 
     531                  araft(ji,jj,jl) = athorn(ji,jj,jl) 
     532               END DO 
     533            END DO 
    421534         END DO 
    422535         ! 
    … …  
    449562      !----------------------------------------------------------------- 
    450563 
    451       aksum(:,:) = athorn(:,:,0) 
     564!$OMP PARALLEL 
     565!$OMP DO schedule(static) private(jj,ji) 
     566      DO jj = 1, jpj  
     567         DO ji = 1, jpi 
     568            aksum(ji,jj) = athorn(ji,jj,0) 
     569         END DO 
     570      END DO 
    452571      ! Transfer function 
    453572      DO jl = 1, jpl !all categories have a specific transfer function 
     573!$OMP DO schedule(static) private(jj,ji,hrmean) 
    454574         DO jj = 1, jpj 
    455575            DO ji = 1, jpi 
    … …  
    476596         END DO 
    477597      END DO 
     598!$OMP END PARALLEL 
    478599      ! 
    479600      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 
     601      CALL wrk_dealloc( jpi,jpj,z_ai ) 
    480602      ! 
    481603   END SUBROUTINE lim_itd_me_ridgeprep 
    … …  
    539661      ! 1) Compute change in open water area due to closing and opening. 
    540662      !------------------------------------------------------------------------------- 
     663!$OMP PARALLEL DO schedule(static) private(jj,ji) 
    541664      DO jj = 1, jpj 
    542665         DO ji = 1, jpi 
    … …  
    568691         END DO 
    569692 
     693!$OMP PARALLEL 
     694!$OMP DO schedule(static) private(ij,jj,ji) 
    570695         DO ij = 1, icells 
    571696            ji = indxi(ij) ; jj = indxj(ij) 
    … …  
    660785         !-------------------------------------------------------------------- 
    661786         DO jk = 1, nlay_i 
     787!$OMP DO schedule(static) private(ij,jj,ji) 
    662788            DO ij = 1, icells 
    663789               ji = indxi(ij) ; jj = indxj(ij) 
    … …  
    687813         DO jl2  = 1, jpl  
    688814            ! over categories to which ridged/rafted ice is transferred 
     815!$OMP DO schedule(static) private(ij,jj,ji,hL,hR,farea) 
    689816            DO ij = 1, icells 
    690817               ji = indxi(ij) ; jj = indxj(ij) 
    … …  
    721848            ! Transfer ice energy to category jl2 by ridging 
    722849            DO jk = 1, nlay_i 
     850!$OMP DO schedule(static) private(ij,jj,ji) 
    723851               DO ij = 1, icells 
    724852                  ji = indxi(ij) ; jj = indxj(ij) 
    … …  
    728856            ! 
    729857         END DO ! jl2 
     858!$OMP END PARALLEL 
    730859          
    731860      END DO ! jl1 (deforming categories) 
    732  
    733861      ! 
    734862      CALL wrk_dealloc( jpij,        indxi, indxj ) 
    … …  
    769897      ! 1) Initialize 
    770898      !------------------------------------------------------------------------------! 
    771       strength(:,:) = 0._wp 
     899!$OMP PARALLEL DO schedule(static) private(jj,ji) 
     900      DO jj = 1, jpj 
     901         DO ji = 1, jpi 
     902            strength(ji,jj) = 0._wp 
     903         END DO 
     904      END DO 
    772905 
    773906      !------------------------------------------------------------------------------! 
    … …  
    781914      IF( kstrngth == 1 ) THEN 
    782915         z1_3 = 1._wp / 3._wp 
     916!$OMP PARALLEL 
    783917         DO jl = 1, jpl 
     918!$OMP DO schedule(static) private(jj,ji) 
    784919            DO jj= 1, jpj 
    785920               DO ji = 1, jpi 
    … …  
    810945         END DO 
    811946    
    812          strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 
     947!$OMP DO schedule(static) private(jj,ji) 
     948         DO jj= 1, jpj 
     949            DO ji = 1, jpi 
     950               strength(ji,jj) = rn_pe_rdg * Cp * strength(ji,jj) / aksum(ji,jj) * tmask(ji,jj,1) 
     951            END DO 
     952         END DO 
     953!$OMP END PARALLEL 
    813954                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 
    814955         ksmooth = 1 
    … …  
    818959      !------------------------------------------------------------------------------! 
    819960      ELSE                      ! kstrngth ne 1:  Hibler (1979) form 
    820          ! 
    821          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1) 
     961!$OMP PARALLEL DO schedule(static) private(jj,ji) 
     962         DO jj= 1, jpj 
     963            DO ji = 1, jpi 
     964               ! 
     965               strength(ji,jj) = rn_pstar * vt_i(ji,jj) * EXP( - rn_crhg * ( 1._wp - at_i(ji,jj) )  ) * tmask(ji,jj,1) 
     966            END DO 
     967         END DO 
    822968         ! 
    823969         ksmooth = 1 
    … …  
    830976      ! CAN BE REMOVED 
    831977      IF( ln_icestr_bvf ) THEN 
     978!$OMP PARALLEL DO schedule(static) private(jj,ji) 
    832979         DO jj = 1, jpj 
    833980            DO ji = 1, jpi 
    … …  
    846993      IF ( ksmooth == 1 ) THEN 
    847994 
     995!$OMP PARALLEL 
     996!$OMP DO schedule(static) private(jj,ji) 
    848997         DO jj = 2, jpjm1 
    849998            DO ji = 2, jpim1 
    … …  
    8591008         END DO 
    8601009 
     1010!$OMP DO schedule(static) private(jj,ji) 
    8611011         DO jj = 2, jpjm1 
    8621012            DO ji = 2, jpim1 
    … …  
    8641014            END DO 
    8651015         END DO 
     1016!$OMP END PARALLEL 
    8661017         CALL lbc_lnk( strength, 'T', 1. ) 
    8671018 
    … …  
    8741025 
    8751026         IF ( numit == nit000 + nn_fsbc - 1 ) THEN 
    876             zstrp1(:,:) = 0._wp 
    877             zstrp2(:,:) = 0._wp 
     1027!$OMP PARALLEL DO schedule(static) private(jj,ji) 
     1028            DO jj = 1, jpj 
     1029               DO ji = 1, jpi 
     1030                  zstrp1(ji,jj) = 0._wp 
     1031                  zstrp2(ji,jj) = 0._wp 
     1032               END DO 
     1033            END DO 
    8781034         ENDIF 
    8791035 
     1036!$OMP PARALLEL DO schedule(static) private(jj,ji,numts_rm,zp) 
    8801037         DO jj = 2, jpjm1 
    8811038            DO ji = 2, jpim1 
Note: See TracChangeset for help on using the changeset viewer.