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 6101 for branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90 – NEMO

Ignore:
Timestamp:
2015-12-17T16:48:41+01:00 (8 years ago)
Author:
cbricaud
Message:

correction of bugs from last update and improvments for CRS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90

    r5602 r6101  
    3838   USE crslbclnk 
    3939   USE lib_mpp 
    40     
     40!cbr   USE ieee_arithmetic    
    4141 
    4242   IMPLICIT NONE 
     
    5757#  include "domzgr_substitute.h90" 
    5858    
    59    !! $Id$ 
    6059CONTAINS 
    6160 
     
    6766      INTEGER  ::  iji, ijj 
    6867      REAL(wp) ::  zmask 
     68      INTEGER  :: ir,jr 
    6969       
    7070      ! Initialize 
     
    130130               !write(narea+5000,*)"mask ",ijie,ijis,ijjs,ijje 
    131131               !ENDIF 
     132 
     133               ir=303-nimpp_crs+1 ; jr=302-njmpp_crs+1 
     134               IF( ji==ir .AND. jj==jr )THEN 
     135                   WRITE(narea+2000,*)"mask",ir,jr,ijis+nimpp-1,ijjs+njmpp-1 
     136               ENDIF 
    132137 
    133138               !IF( ijje .GT. jpj )WRITE(narea+200,*)"BUG",jj,ijjs,ijje,SHAPE(tmask) ; call flush(narea+200) 
     
    215220      INTEGER :: ji, jj, jk                   ! dummy loop indices 
    216221      INTEGER :: ijis, ijjs 
     222      INTEGER  :: ir,jr 
    217223 
    218224   
     
    225231                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs) 
    226232                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs) 
     233                  ir=303-nimpp_crs+1 ; jr=302-njmpp_crs+1 
     234                  WRITE(narea+2000,*)"coordT1",ir,jr 
     235                  IF( ji==ir .AND. jj==jr )THEN 
     236                     WRITE(narea+2000,*)"coordT",ir,jr,ijis+nimpp-1,ijjs+njmpp-1 
     237                  ENDIF 
    227238               ENDDO 
    228239            ENDDO 
     
    530541      !!  
    531542      !!  Arguments 
    532       REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_fld   ! T, U, V or W on parent grid 
    533       CHARACTER(len=3),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN 
     543      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)        :: p_fld   ! T, U, V or W on parent grid 
     544      CHARACTER(len=*),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN 
    534545      CHARACTER(len=1),                         INTENT(in)           :: cd_type    ! grid type U,V  
    535546      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_mask    ! Parent grid T,U,V mask 
     
    547558      INTEGER  :: ii, ij, ijie, ijje, je_2 
    548559      REAL(wp) :: zflcrs, zsfcrs    
    549       REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf, zsurfmsk, zmask   
     560      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf, zsurfmsk, zmask,ztabtmp 
    550561      INTEGER  :: iji, ijj 
     562      INTEGER  :: ir,jr 
     563      REAL(wp), DIMENSION(nn_factx,nn_facty):: ztmp 
     564      REAL(wp), DIMENSION(nn_factx*nn_facty):: ztmp1 
     565      REAL(wp), DIMENSION(:), ALLOCATABLE   :: ztmp2 
     566      INTEGER , DIMENSION(1)  :: zdim1 
     567      REAL(wp) :: zmin,zmax 
    551568      !!----------------------------------------------------------------   
    552569    
     
    554571 
    555572      SELECT CASE ( cd_op ) 
    556        
     573  
    557574         CASE ( 'VOL' ) 
    558575       
     
    633650                             &     + p_fld(ji+1,jj+2,jk) * zsurfmsk(ji+1,jj+2,jk) & 
    634651                             &     + p_fld(ji+2,jj+2,jk) * zsurfmsk(ji+2,jj+2,jk)  
    635  
    636652                           zsfcrs =  zsurf(ji,jj  ,jk) + zsurf(ji+1,jj  ,jk) + zsurf(ji+2,jj  ,jk) & 
    637653                             &     + zsurf(ji,jj+1,jk) + zsurf(ji+1,jj+1,jk) + zsurf(ji+2,jj+1,jk) & 
    638654                             &     + zsurf(ji,jj+2,jk) + zsurf(ji+1,jj+2,jk) + zsurf(ji+2,jj+2,jk)  
    639655                            ! 
     656!cbr                            IF( ieee_is_nan(p_fld_crs(ii,ij,jk))) THEN 
     657 
    640658                           p_fld_crs(ii,ij,jk) = zflcrs 
    641659                           IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,jk) = zflcrs / zsfcrs 
     
    648666 
    649667              CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk ) 
    650  
     668         CASE ( 'LOGVOL' ) 
     669 
     670            CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk, ztabtmp ) 
     671 
     672            zmin=MINVAL(p_fld) ; zmax=MAXVAL(p_fld);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"p_fld",zmin,zmax; CALL flush(numout) 
     673 
     674            ztabtmp(:,:,:)=0._wp 
     675            WHERE(p_fld* p_mask .NE. 0._wp ) ztabtmp =  LOG10(p_fld * p_mask)*p_mask 
     676            zmin=MINVAL(ztabtmp) ; zmax=MAXVAL(ztabtmp);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"LOG()",zmin,zmax; CALL flush(numout) 
     677            ztabtmp = ztabtmp * p_mask 
     678            zmin=MINVAL(ztabtmp) ; zmax=MAXVAL(ztabtmp);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"LOG()*tmask",zmin,zmax; CALL flush(numout) 
     679 
     680            SELECT CASE ( cd_type ) 
     681 
     682               CASE( 'T', 'W' ) 
     683                     DO jk = 1, jpk 
     684                        zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk) 
     685                        zsurfmsk(:,:,jk) =  zsurf(:,:,jk) 
     686                    ENDDO 
     687 
     688                  IF( nldj_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2 
     689                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN 
     690                        je_2 = mje_crs(2) 
     691                        DO jk = 1, jpk 
     692                           DO ji = nistr, niend, nn_factx 
     693                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2 
     694                              zflcrs =  ztabtmp(ji  ,je_2,jk) * zsurfmsk(ji  ,je_2,jk)   & 
     695                                &     + ztabtmp(ji+1,je_2,jk) * zsurfmsk(ji+1,je_2,jk)   & 
     696                                &     + ztabtmp(ji+2,je_2,jk) * zsurfmsk(ji+2,je_2,jk) 
     697 
     698                              zsfcrs =  zsurf(ji,je_2,jk) + zsurf(ji+1,je_2,jk) + zsurf(ji+2,je_2,jk) 
     699                              ! 
     700                              p_fld_crs(ii,2,jk) = 0._wp 
     701                              IF( zsfcrs /= 0.0 )  p_fld_crs(ii,2,jk) = zflcrs / zsfcrs 
     702                              p_fld_crs(ii,2,jk) = 10 ** ( p_fld_crs(ii,2,jk) * p_mask_crs(ii,2,jk) ) * p_mask_crs(ii,2,jk) 
     703                           ENDDO 
     704                        ENDDO 
     705                     ENDIF 
     706                  ELSE 
     707                     je_2 = mjs_crs(2) 
     708                     DO jk = 1, jpk 
     709                        DO ji = nistr, niend, nn_factx 
     710                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2 
     711                           zflcrs =  ztabtmp(ji  ,je_2  ,jk) * zsurfmsk(ji  ,je_2  ,jk) & 
     712                             &     + ztabtmp(ji+1,je_2  ,jk) * zsurfmsk(ji+1,je_2  ,jk) & 
     713                             &     + ztabtmp(ji+2,je_2  ,jk) * zsurfmsk(ji+2,je_2  ,jk) & 
     714                             &     + ztabtmp(ji  ,je_2+1,jk) * zsurfmsk(ji  ,je_2+1,jk) & 
     715                             &     + ztabtmp(ji+1,je_2+1,jk) * zsurfmsk(ji+1,je_2+1,jk) & 
     716                             &     + ztabtmp(ji+2,je_2+1,jk) * zsurfmsk(ji+2,je_2+1,jk) & 
     717                             &     + ztabtmp(ji  ,je_2+2,jk) * zsurfmsk(ji  ,je_2+2,jk) & 
     718                             &     + ztabtmp(ji+1,je_2+2,jk) * zsurfmsk(ji+1,je_2+2,jk) & 
     719                             &     + ztabtmp(ji+2,je_2+2,jk) * zsurfmsk(ji+2,je_2+2,jk) 
     720 
     721                           zsfcrs =  zsurf(ji,je_2  ,jk) + zsurf(ji+1,je_2  ,jk) + zsurf(ji+2,je_2  ,jk) & 
     722                             &     + zsurf(ji,je_2+1,jk) + zsurf(ji+1,je_2+1,jk) + zsurf(ji+2,je_2+1,jk) & 
     723                             &     + zsurf(ji,je_2+2,jk) + zsurf(ji+1,je_2+2,jk) + zsurf(ji+2,je_2+2,jk) 
     724                            ! 
     725                            p_fld_crs(ii,2,jk) = 0._wp 
     726                            IF( zsfcrs /= 0.0 )  p_fld_crs(ii,2,jk) = zflcrs / zsfcrs 
     727                            p_fld_crs(ii,2,jk) = 10 ** ( p_fld_crs(ii,2,jk) * p_mask_crs(ii,2,jk) ) * p_mask_crs(ii,2,jk) 
     728                        ENDDO 
     729                     ENDDO 
     730                  ENDIF 
     731                  ! 
     732                  DO jk = 1, jpk 
     733                     DO jj  = njstr, njend, nn_facty 
     734                        DO ji = nistr, niend, nn_factx 
     735                           ii = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid 
     736                           ij = ( jj - njstr ) * rfacty_r + 3 
     737                           zflcrs =  ztabtmp(ji  ,jj  ,jk) * zsurfmsk(ji  ,jj  ,jk) & 
     738                             &     + ztabtmp(ji+1,jj  ,jk) * zsurfmsk(ji+1,jj  ,jk) & 
     739                             &     + ztabtmp(ji+2,jj  ,jk) * zsurfmsk(ji+2,jj  ,jk) & 
     740                             &     + ztabtmp(ji  ,jj+1,jk) * zsurfmsk(ji  ,jj+1,jk) & 
     741                             &     + ztabtmp(ji+1,jj+1,jk) * zsurfmsk(ji+1,jj+1,jk) & 
     742                             &     + ztabtmp(ji+2,jj+1,jk) * zsurfmsk(ji+2,jj+1,jk) & 
     743                             &     + ztabtmp(ji  ,jj+2,jk) * zsurfmsk(ji  ,jj+2,jk) & 
     744                             &     + ztabtmp(ji+1,jj+2,jk) * zsurfmsk(ji+1,jj+2,jk) & 
     745                             &     + ztabtmp(ji+2,jj+2,jk) * zsurfmsk(ji+2,jj+2,jk) 
     746                           zsfcrs =  zsurf(ji,jj  ,jk) + zsurf(ji+1,jj  ,jk) + zsurf(ji+2,jj  ,jk) & 
     747                             &     + zsurf(ji,jj+1,jk) + zsurf(ji+1,jj+1,jk) + zsurf(ji+2,jj+1,jk) & 
     748                             &     + zsurf(ji,jj+2,jk) + zsurf(ji+1,jj+2,jk) + zsurf(ji+2,jj+2,jk) 
     749                            ! 
     750                           p_fld_crs(ii,ij,jk) = 0._wp 
     751                           IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,jk) = zflcrs / zsfcrs 
     752                           p_fld_crs(ii,ij,jk) = 10 ** ( p_fld_crs(ii,ij,jk) *  p_mask_crs(ii,ij,jk) ) * p_mask_crs(ii,ij,jk) 
     753                        ENDDO 
     754                     ENDDO 
     755                  ENDDO 
     756               CASE DEFAULT 
     757                    STOP 
     758               END SELECT 
     759 
     760 
     761              !WHERE( p_fld .NE. 0._wp ) p_fld=10**(p_fld) 
     762              !zmin=MINVAL(p_fld) ; zmax=MAXVAL(p_fld);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld)",zmin,zmax ; CALL flush(numout) 
     763              !p_fld = p_fld * p_mask 
     764              !zmin=MINVAL(p_fld) ; zmax=MAXVAL(p_fld);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld)*tmask",zmin,zmax ; CALL flush(numout) 
     765 
     766              zmin=MINVAL(p_fld_crs) ; zmax=MAXVAL(p_fld_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"p_fld_crs",zmin,zmax; CALL flush(numout) 
     767              !p_fld_crs=10**(p_fld_crs*p_mask_crs) 
     768              !zmin=MINVAL(p_fld_crs) ; zmax=MAXVAL(p_fld_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld_crs)",zmin,zmax; CALL flush(numout) 
     769              !p_fld_crs=p_fld_crs*p_mask_crs 
     770              !zmin=MINVAL(p_fld_crs) ; zmax=MAXVAL(p_fld_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld_crs)*tmask",zmin,zmax; CALL flush(numout) 
     771 
     772              CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk ,ztabtmp ) 
     773         CASE ( 'MED' ) 
     774 
     775            CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk ) 
     776 
     777            SELECT CASE ( cd_type ) 
     778 
     779               CASE( 'T', 'W' ) 
     780                     DO jk = 1, jpk 
     781                        zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk) 
     782                        zsurfmsk(:,:,jk) =  zsurf(:,:,jk) 
     783                    ENDDO 
     784 
     785                  IF( nldj_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2 
     786                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN 
     787                        je_2 = mje_crs(2) 
     788                        DO jk = 1, jpk 
     789                           DO ji = nistr, niend, nn_factx 
     790                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2 
     791 
     792                              ztmp1(:) = 0._wp 
     793                              ztmp1(1:3) =  p_fld(ji:ji+2,je_2,jk) 
     794                              CALL PIKSRT(nn_factx*nn_facty,ztmp1) 
     795                              ir=0 
     796                              jr=1 
     797                              DO WHILE( jr .LE. nn_factx*nn_facty ) 
     798                                 IF( ztmp1(jr) == 0. )THEN 
     799                                    ir=jr 
     800                                    jr=jr+1 
     801                                 ELSE 
     802                                    EXIT 
     803                                 ENDIF 
     804                              ENDDO 
     805                              IF( ir .LE. nn_factx*nn_facty-1 )THEN 
     806                                 ALLOCATE( ztmp2(nn_factx*nn_facty-ir) ) 
     807                                 ztmp2(1:nn_factx*nn_facty-ir) = ztmp1(1+ir:nn_factx*nn_facty) 
     808                                 jr=INT( 0.5 * REAL(nn_factx*nn_facty-ir,wp) )+1 
     809                                 p_fld_crs(ii,2,jk) = ztmp2(jr) 
     810                                 DEALLOCATE( ztmp2 ) 
     811                              ELSE 
     812                                 p_fld_crs(ii,ij,jk) = 0._wp 
     813                              ENDIF 
     814 
     815                           ENDDO 
     816                        ENDDO 
     817                     ENDIF 
     818                  ELSE 
     819                     je_2 = mjs_crs(2) 
     820                     DO jk = 1, jpk 
     821                        DO ji = nistr, niend, nn_factx 
     822                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2 
     823                            
     824                           ztmp(:,:)= p_fld(ji:ji+2,je_2:je_2+2,jk) 
     825                           zdim1(1)=nn_factx*nn_facty 
     826                           ztmp1(:) = RESHAPE( ztmp(:,:) , zdim1 ) 
     827                           CALL PIKSRT(nn_factx*nn_facty,ztmp1) 
     828                           ir=0 
     829                           jr=1 
     830                           DO WHILE( jr .LE. nn_factx*nn_facty ) 
     831                              IF( ztmp1(jr) == 0. ) THEN 
     832                                 ir=jr 
     833                                 jr=jr+1 
     834                              ELSE 
     835                                 EXIT 
     836                              ENDIF 
     837                           ENDDO 
     838                           IF( ir .LE. nn_factx*nn_facty-1 )THEN 
     839                              ALLOCATE( ztmp2(nn_factx*nn_facty-ir) ) 
     840                              ztmp2(1:nn_factx*nn_facty-ir) = ztmp1(1+ir:nn_factx*nn_facty) 
     841                              jr=INT( 0.5 * REAL(nn_factx*nn_facty-ir,wp) )+1 
     842                              p_fld_crs(ii,2,jk) = ztmp2(jr) 
     843                              DEALLOCATE( ztmp2 ) 
     844                           ELSE 
     845                           p_fld_crs(ii,ij,jk) = 0._wp 
     846                           ENDIF 
     847 
     848                        ENDDO 
     849                     ENDDO 
     850                  ENDIF 
     851                  ! 
     852                  DO jk = 1, jpk 
     853                     DO jj  = njstr, njend, nn_facty 
     854                        DO ji = nistr, niend, nn_factx 
     855                           ii = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid 
     856                           ij = ( jj - njstr ) * rfacty_r + 3 
     857 
     858                           ztmp(:,:)= p_fld(ji:ji+2,jj:jj+2,jk)  
     859                           zdim1(1)=nn_factx*nn_facty 
     860                           ztmp1(:) = RESHAPE( ztmp(:,:) , zdim1 ) 
     861                           CALL PIKSRT(nn_factx*nn_facty,ztmp1) 
     862                           ir=0 
     863                           jr=1 
     864                           DO WHILE( jr .LE. nn_factx*nn_facty ) 
     865                              IF( ztmp1(jr) == 0. ) THEN 
     866                                 ir=jr 
     867                                 jr=jr+1 
     868                              ELSE 
     869                                 EXIT 
     870                              ENDIF 
     871                           ENDDO 
     872                           IF( ir .LE. nn_factx*nn_facty-1 )THEN 
     873                              ALLOCATE( ztmp2(nn_factx*nn_facty-ir) ) 
     874                              ztmp2(1:nn_factx*nn_facty-ir) = ztmp1(1+ir:nn_factx*nn_facty) 
     875                              jr=INT( 0.5 * REAL(nn_factx*nn_facty-ir,wp) )+1 
     876                              p_fld_crs(ii,ij,jk) = ztmp2(jr) 
     877                              DEALLOCATE( ztmp2 ) 
     878                           ELSE 
     879                              p_fld_crs(ii,ij,jk) = 0._wp 
     880                           ENDIF 
     881 
     882                        ENDDO 
     883                     ENDDO 
     884                  ENDDO 
     885               CASE DEFAULT 
     886                    STOP 
     887               END SELECT 
     888 
     889              CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk ) 
     890  
    651891         CASE ( 'SUM' ) 
    652892          
     
    23902630      CALL crs_lbc_lnk( p_surf_crs_msk, cd_type, 1.0, pval=1.0 ) 
    23912631 
    2392       CALL wrk_dealloc( jpi, jpj, jpk, zsurfmsk, zsurf ) 
     2632      CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk ) 
    23932633 
    23942634   END SUBROUTINE crs_dom_sfc 
     
    28933133      ENDDO 
    28943134      
     3135      CALL wrk_alloc( jpi_crs, jpj_crs, zmbk ) 
     3136 
    28953137      zmbk(:,:) = 0.0 
    28963138      zmbk(:,:) = REAL( mbathy_crs(:,:), wp ) ;   CALL crs_lbc_lnk(zmbk,'T',1.0)   ;   mbathy_crs(:,:) = INT( zmbk(:,:) ) 
     
    29213163   END SUBROUTINE crs_dom_bat 
    29223164 
     3165   SUBROUTINE PIKSRT(N,ARR) 
     3166     INTEGER                  ,INTENT(IN) :: N 
     3167     REAL(kind=8),DIMENSION(N),INTENT(INOUT) :: ARR 
     3168 
     3169     INTEGER      :: i,j 
     3170     REAL(kind=8) :: a 
     3171     !!---------------------------------------------------------------- 
     3172 
     3173     DO j=2, N 
     3174       a=ARR(j) 
     3175       DO i=j-1,1,-1 
     3176          IF(ARR(i)<=a) goto 10 
     3177          ARR(i+1)=ARR(i) 
     3178       ENDDO 
     3179       i=0 
     318010     ARR(i+1)=a 
     3181     ENDDO 
     3182     RETURN 
     3183 
     3184   END SUBROUTINE PIKSRT 
     3185 
    29233186 
    29243187END MODULE crsdom 
Note: See TracChangeset for help on using the changeset viewer.