Changeset 845


Ignore:
Timestamp:
05/04/19 21:35:04 (5 years ago)
Author:
dubos
Message:

devel : reorganized dissip_gcm.f90

File:
1 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/src/dissip/dissip_gcm.f90

    r714 r845  
    11MODULE dissip_gcm_mod 
    22  USE icosa 
    3  
     3  USE omp_para 
     4  USE trace 
     5  IMPLICIT NONE 
    46  PRIVATE 
    57 
    … …  
    4143  REAL,SAVE    :: dtdissip 
    4244!$OMP THREADPRIVATE(dtdissip) 
    43   
     45   
    4446  PUBLIC init_dissip, dissip 
     47 
    4548CONTAINS 
    46  
     49   
    4750  SUBROUTINE allocate_dissip 
    48   USE icosa 
    49   IMPLICIT NONE   
    5051    CALL allocate_field(f_due_diss1,field_u,type_real,llm) 
    5152    CALL allocate_field(f_due_diss2,field_u,type_real,llm) 
    … …  
    5859  
    5960  SUBROUTINE init_dissip 
    60   USE icosa 
    61   USE disvert_mod 
    62   USE mpi_mod 
    63   USE mpipara 
    64   USE transfert_mod 
    65   USE time_mod 
    66   USE transfert_omp_mod 
    67   USE omp_para 
    68   IMPLICIT NONE 
    69    
    70   TYPE(t_field),POINTER,SAVE  :: f_u(:) 
    71   TYPE(t_field),POINTER,SAVE  :: f_du(:) 
    72   REAL(rstd),POINTER     :: u(:) 
    73   REAL(rstd),POINTER     :: du(:) 
    74   TYPE(t_field),POINTER,SAVE  :: f_theta(:) 
    75   TYPE(t_field),POINTER ,SAVE :: f_dtheta(:) 
    76   REAL(rstd),POINTER    :: theta(:) 
    77   REAL(rstd),POINTER    :: dtheta(:) 
    78   REAL(rstd)            :: dumax,dumax1 
    79   REAL(rstd)            :: dthetamax,dthetamax1 
    80   REAL(rstd)            :: r 
    81   REAL(rstd)            :: tau 
    82   REAL(rstd)            :: zz, zvert, fact 
    83   INTEGER               :: l 
    84   CHARACTER(len=255)    :: rayleigh_friction_key 
    85   REAL(rstd)            :: mintau 
    86   INTEGER               :: seed_size 
    87   INTEGER,ALLOCATABLE   :: seed(:) 
    88    
    89              
    90   INTEGER :: i,j,ij,ind,it,iter,M 
    91  
    92    rayleigh_friction_key='none' 
    93    CALL getin("rayleigh_friction_type",rayleigh_friction_key) 
    94    SELECT CASE(TRIM(rayleigh_friction_key)) 
    95    CASE('none') 
    96       rayleigh_friction_type=0 
    97       IF (is_master) PRINT *, 'No Rayleigh friction' 
    98    CASE('dcmip2_schaer_noshear') 
    99       rayleigh_friction_type=1 
    100       rayleigh_shear=0 
    101       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1' 
    102    CASE('dcmip2_schaer_shear') 
    103       rayleigh_shear=1 
    104       rayleigh_friction_type=2 
    105       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2' 
    106    CASE('giant_liu_schneider')  
    107       rayleigh_friction_type=99 
    108       IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010' 
    109    CASE DEFAULT 
    110       IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 
    111         ' in dissip_gcm.f90/init_dissip' 
    112       STOP 
    113    END SELECT 
    114  
    115    IF(rayleigh_friction_type>0) THEN 
    116       rayleigh_tau=0. 
    117       CALL getin("rayleigh_friction_tau",rayleigh_tau) 
    118       rayleigh_tau = rayleigh_tau / scale_factor 
    119       IF(rayleigh_tau<=0) THEN 
    120          IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau 
    121          STOP 
    122       END IF 
    123       IF(rayleigh_friction_type == 99) THEN 
    124          rayleigh_limlat=0. 
    125          CALL getin("rayleigh_limlat",rayleigh_limlat) 
    126          rayleigh_limlat = rayleigh_limlat*3.14159/180. 
    127       ENDIF 
    128    END IF 
    129  
     61    REAL(rstd)            :: tau   
     62    CHARACTER(len=255)    :: rayleigh_friction_key 
     63    rayleigh_friction_key='none' 
     64    CALL getin("rayleigh_friction_type",rayleigh_friction_key) 
     65    SELECT CASE(TRIM(rayleigh_friction_key)) 
     66    CASE('none') 
     67       rayleigh_friction_type=0 
     68       IF (is_master) PRINT *, 'No Rayleigh friction' 
     69    CASE('dcmip2_schaer_noshear') 
     70       rayleigh_friction_type=1 
     71       rayleigh_shear=0 
     72       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1' 
     73    CASE('dcmip2_schaer_shear') 
     74       rayleigh_shear=1 
     75       rayleigh_friction_type=2 
     76       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2' 
     77    CASE('giant_liu_schneider')  
     78       rayleigh_friction_type=99 
     79       IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010' 
     80    CASE DEFAULT 
     81       IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 
     82            ' in dissip_gcm.f90/init_dissip' 
     83       STOP 
     84    END SELECT 
     85     
     86    IF(rayleigh_friction_type>0) THEN 
     87       rayleigh_tau=0. 
     88       CALL getin("rayleigh_friction_tau",rayleigh_tau) 
     89       rayleigh_tau = rayleigh_tau / scale_factor 
     90       IF(rayleigh_tau<=0) THEN 
     91          IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau 
     92          STOP 
     93       END IF 
     94       IF(rayleigh_friction_type == 99) THEN 
     95          rayleigh_limlat=0. 
     96          CALL getin("rayleigh_limlat",rayleigh_limlat) 
     97          rayleigh_limlat = rayleigh_limlat*3.14159/180. 
     98       ENDIF 
     99    END IF 
     100     
    130101    CALL allocate_dissip 
     102 
     103    CALL init_message(f_due_diss1,req_e1_vect,req_due) 
     104    CALL init_message(f_dtheta_diss,req_i1,req_dtheta) 
     105 
     106    tau_graddiv(:)=5000 
     107    CALL getin("tau_graddiv",tau) 
     108    tau_graddiv(:)=tau/scale_factor 
     109 
     110    CALL getin("nitergdiv",nitergdiv) 
     111     
     112    tau_gradrot(:)=5000 
     113    CALL getin("tau_gradrot",tau) 
     114    tau_gradrot(:)=tau/scale_factor 
     115 
     116    CALL getin("nitergrot",nitergrot) 
     117     
     118 
     119    tau_divgrad(:)=5000 
     120    CALL getin("tau_divgrad",tau) 
     121    tau_divgrad(:)=tau/scale_factor 
     122 
     123    CALL getin("niterdivgrad",niterdivgrad) 
     124 
     125    CALL dissip_constants 
     126    CALL dissip_profile 
     127  END SUBROUTINE init_dissip 
     128 
     129  SUBROUTINE dissip_constants 
     130    ! The SAVE attribute is required here, otherwise allocate_field will not work with OpenMP 
     131    TYPE(t_field),POINTER, SAVE  :: f_u(:) 
     132    TYPE(t_field),POINTER, SAVE  :: f_du(:) 
     133    TYPE(t_field),POINTER, SAVE  :: f_theta(:) 
     134    TYPE(t_field),POINTER, SAVE  :: f_dtheta(:) 
     135    REAL(rstd),POINTER :: u(:) 
     136    REAL(rstd),POINTER :: du(:) 
     137    REAL(rstd),POINTER :: theta(:) 
     138    REAL(rstd),POINTER :: dtheta(:) 
     139    REAL(rstd) :: dumax, dthetamax 
     140    INTEGER :: it, iter, ind 
    131141    CALL allocate_field(f_u,field_u,type_real) 
    132142    CALL allocate_field(f_du,field_u,type_real) 
    133143    CALL allocate_field(f_theta,field_t,type_real) 
    134144    CALL allocate_field(f_dtheta,field_t,type_real) 
    135      
    136     CALL init_message(f_due_diss1,req_e1_vect,req_due) 
    137     CALL init_message(f_dtheta_diss,req_i1,req_dtheta) 
    138  
    139     tau_graddiv(:)=5000 
    140     CALL getin("tau_graddiv",tau) 
    141     tau_graddiv(:)=tau/scale_factor 
    142  
    143     CALL getin("nitergdiv",nitergdiv) 
    144      
    145     tau_gradrot(:)=5000 
    146     CALL getin("tau_gradrot",tau) 
    147     tau_gradrot(:)=tau/scale_factor 
    148  
    149     CALL getin("nitergrot",nitergrot) 
    150      
    151  
    152     tau_divgrad(:)=5000 
    153     CALL getin("tau_divgrad",tau) 
    154     tau_divgrad(:)=tau/scale_factor 
    155  
    156     CALL getin("niterdivgrad",niterdivgrad) 
    157  
    158  
    159     cgraddiv=-1 
    160     cdivgrad=-1 
    161     cgradrot=-1 
    162  
    163 !$OMP BARRIER 
    164 !$OMP MASTER 
    165     DO ind=1,ndomain 
    166       CALL swap_dimensions(ind) 
    167       CALL swap_geometry(ind) 
    168       u=f_u(ind) 
    169  
    170 ! set random seed to get reproductibility when using a different number of process 
    171       CALL RANDOM_SEED(size=M) 
    172       CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
    173  
    174       DO j=jj_begin,jj_end 
    175         DO i=ii_begin,ii_end 
    176           ij=(j-1)*iim+i    
    177           CALL RANDOM_NUMBER(r) 
    178           u(ij+u_right)=r-0.5 
    179           CALL RANDOM_NUMBER(r) 
    180           u(ij+u_lup)=r-0.5 
    181           CALL RANDOM_NUMBER(r) 
    182           u(ij+u_ldown)=r-0.5 
    183         ENDDO 
    184       ENDDO         
    185     ENDDO 
    186 !$OMP END MASTER   
    187 !$OMP BARRIER 
    188  
     145 
     146    PRINT *, '=========', SHAPE(f_u) 
     147 
     148    !--------------------------- Compute cgraddiv ---------------------------- 
     149    cgraddiv=-1. 
     150    CALL random_vel(f_u)     
    189151    DO it=1,20 
    190        
    191       dumax=0 
    192       DO iter=1,nitergdiv 
    193         CALL transfert_request(f_u,req_e1_vect) 
    194         DO ind=1,ndomain 
    195           IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    196           CALL swap_dimensions(ind) 
    197           CALL swap_geometry(ind) 
    198           u=f_u(ind) 
    199           du=f_du(ind) 
    200           CALL compute_gradiv(u,du,1,1) 
    201           u=du 
    202         ENDDO 
    203       ENDDO 
    204  
    205       CALL transfert_request(f_du,req_e1_vect) 
    206        
    207       DO ind=1,ndomain 
    208         IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    209         CALL swap_dimensions(ind) 
    210         CALL swap_geometry(ind) 
    211         u=f_u(ind) 
    212         du=f_du(ind) 
    213           
    214         DO j=jj_begin,jj_end 
    215           DO i=ii_begin,ii_end 
    216             ij=(j-1)*iim+i    
    217             if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 
    218             if (le(ij+u_lup)>1e-100)   dumax=MAX(dumax,ABS(du(ij+u_lup))) 
    219             if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 
     152       DO iter=1,nitergdiv 
     153          CALL transfert_request(f_u,req_e1_vect) 
     154          DO ind=1,ndomain 
     155             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     156             CALL swap_dimensions(ind) 
     157             CALL swap_geometry(ind) 
     158             u=f_u(ind) 
     159             du=f_du(ind) 
     160             CALL compute_gradiv(u,du,1,1) 
    220161          ENDDO 
    221         ENDDO 
    222       ENDDO 
    223  
    224       IF (using_mpi) THEN  
    225         CALL reduce_max_omp(dumax,dumax1) 
    226 !$OMP MASTER         
    227         CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    228 !$OMP END MASTER      
    229         CALL bcast_omp(dumax)   
    230       ELSE 
    231         CALL allreduce_max_omp(dumax,dumax1) 
    232         dumax=dumax1 
    233       ENDIF   
    234                          
    235       DO ind=1,ndomain 
    236         IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    237         CALL swap_dimensions(ind) 
    238         CALL swap_geometry(ind) 
    239         u=f_u(ind) 
    240         du=f_du(ind) 
    241         u=du/dumax 
    242       ENDDO 
    243       IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax 
    244  
    245     ENDDO  
    246     IF (is_master) PRINT *,"gradiv : dumax",dumax 
    247     IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., & 
    248                               'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000 
    249     IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', & 
    250                                3./340.*dumax**(-.5/nitergdiv) 
     162       ENDDO 
     163       CALL transfert_request(f_du,req_e1_vect)        
     164       dumax = max_vel(f_du)        
     165       CALL rescale_field(1./dumax, f_du, f_u) 
     166       IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax 
     167    ENDDO 
    251168 
    252169    cgraddiv=dumax**(-1./nitergdiv) 
    253     IF (is_master) PRINT *,"cgraddiv : ",cgraddiv 
    254  
    255 !$OMP BARRIER 
    256 !$OMP MASTER 
    257     DO ind=1,ndomain 
    258       CALL swap_dimensions(ind) 
    259       CALL swap_geometry(ind) 
    260       u=f_u(ind) 
    261  
    262 ! set random seed to get reproductibility when using a different number of process 
    263       CALL RANDOM_SEED(size=M) 
    264       CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
    265   
    266        DO j=jj_begin,jj_end 
    267         DO i=ii_begin,ii_end 
    268           ij=(j-1)*iim+i    
    269           CALL RANDOM_NUMBER(r) 
    270           u(ij+u_right)=r-0.5 
    271           CALL RANDOM_NUMBER(r) 
    272           u(ij+u_lup)=r-0.5 
    273           CALL RANDOM_NUMBER(r) 
    274           u(ij+u_ldown)=r-0.5 
    275         ENDDO 
    276       ENDDO         
    277     ENDDO 
    278 !$OMP END MASTER   
    279 !$OMP BARRIER 
    280  
    281  
     170    IF (is_master) THEN 
     171       PRINT *,"gradiv : dumax",dumax 
     172       PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., & 
     173         'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000 
     174       PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', & 
     175         3./340.*dumax**(-.5/nitergdiv) 
     176       PRINT *,"cgraddiv : ",cgraddiv 
     177    END IF 
     178     
     179    !----------------- Compute cgradrot -------------------- 
     180    cgradrot=-1. 
     181    CALL random_vel(f_u) 
    282182    DO it=1,20 
    283   
    284       dumax=0 
    285       DO iter=1,nitergrot 
    286         CALL transfert_request(f_u,req_e1_vect) 
    287         DO ind=1,ndomain 
    288           IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    289           CALL swap_dimensions(ind) 
    290           CALL swap_geometry(ind) 
    291           u=f_u(ind) 
    292           du=f_du(ind) 
    293           CALL compute_gradrot(u,du,1,1) 
    294           u=du 
    295         ENDDO 
    296       ENDDO 
    297  
    298       CALL transfert_request(f_du,req_e1_vect) 
    299        
    300       DO ind=1,ndomain 
    301         IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    302         CALL swap_dimensions(ind) 
    303         CALL swap_geometry(ind) 
    304         u=f_u(ind) 
    305         du=f_du(ind) 
    306          
    307         DO j=jj_begin,jj_end 
    308           DO i=ii_begin,ii_end 
    309             ij=(j-1)*iim+i    
    310             if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 
    311             if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 
    312             if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 
     183       DO iter=1,nitergrot 
     184          CALL transfert_request(f_u,req_e1_vect) 
     185          DO ind=1,ndomain 
     186             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     187             CALL swap_dimensions(ind) 
     188             CALL swap_geometry(ind) 
     189             u=f_u(ind) 
     190             du=f_du(ind) 
     191             CALL compute_gradrot(u,du,1,1) 
     192             u=du 
    313193          ENDDO 
    314         ENDDO 
    315       ENDDO 
    316  
    317       IF (using_mpi) THEN  
    318         CALL reduce_max_omp(dumax,dumax1) 
    319 !$OMP MASTER         
    320         CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    321 !$OMP END MASTER      
    322         CALL bcast_omp(dumax)   
    323       ELSE 
    324         CALL allreduce_max_omp(dumax,dumax1) 
    325         dumax=dumax1 
    326       ENDIF   
    327  
    328        
    329       DO ind=1,ndomain 
    330         IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    331         CALL swap_dimensions(ind) 
    332         CALL swap_geometry(ind) 
    333         u=f_u(ind) 
    334         du=f_du(ind) 
    335         u=du/dumax 
    336       ENDDO 
    337        
    338       IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax 
    339  
    340     ENDDO  
    341     IF (is_master) PRINT *,"gradrot : dumax",dumax 
    342    
     194       ENDDO 
     195       CALL transfert_request(f_du,req_e1_vect) 
     196       dumax = max_vel(f_du) 
     197       CALL rescale_field(1./dumax, f_du, f_u) 
     198       IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax 
     199    ENDDO 
     200 
    343201    cgradrot=dumax**(-1./nitergrot)  
     202    IF (is_master) PRINT *,"gradrot : dumax",dumax   
    344203    IF (is_master) PRINT *,"cgradrot : ",cgradrot 
    345204    
    346  
    347  
    348 !$OMP BARRIER 
    349 !$OMP MASTER 
    350     DO ind=1,ndomain 
    351       CALL swap_dimensions(ind) 
    352       CALL swap_geometry(ind) 
    353       theta=f_theta(ind) 
    354  
    355 ! set random seed to get reproductibility when using a different number of process 
    356       CALL RANDOM_SEED(size=M) 
    357       CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
    358  
    359       DO j=jj_begin,jj_end 
    360         DO i=ii_begin,ii_end 
    361           ij=(j-1)*iim+i    
    362           CALL RANDOM_NUMBER(r) 
    363           theta(ij)=r-0.5 
    364         ENDDO 
    365       ENDDO   
    366     ENDDO 
    367 !$OMP END MASTER 
    368 !$OMP BARRIER 
    369  
    370     DO it=1,20 
    371  
    372       dthetamax=0 
    373       DO iter=1,niterdivgrad 
    374         CALL transfert_request(f_theta,req_i1) 
    375         DO ind=1,ndomain 
    376           IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    377           CALL swap_dimensions(ind) 
    378           CALL swap_geometry(ind) 
    379           theta=f_theta(ind) 
    380           dtheta=f_dtheta(ind) 
    381           CALL compute_divgrad(theta,dtheta,1,1) 
    382           theta=dtheta 
    383         ENDDO 
    384       ENDDO 
    385  
    386       CALL transfert_request(f_dtheta,req_i1) 
    387        
    388       DO ind=1,ndomain 
    389         IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    390         CALL swap_dimensions(ind) 
    391         CALL swap_geometry(ind) 
    392         theta=f_theta(ind) 
    393         dtheta=f_dtheta(ind) 
    394          
    395         DO j=jj_begin,jj_end 
    396           DO i=ii_begin,ii_end 
    397             ij=(j-1)*iim+i    
    398             dthetamax=MAX(dthetamax,ABS(dtheta(ij))) 
     205    !----------------- Compute cdivgrad -------------------- 
     206 
     207    cdivgrad=-1. 
     208    CALL random_scalar(f_theta)     
     209    DO it=1,20        
     210       DO iter=1,niterdivgrad 
     211          CALL transfert_request(f_theta,req_i1) 
     212          DO ind=1,ndomain 
     213             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     214             CALL swap_dimensions(ind) 
     215             CALL swap_geometry(ind) 
     216             theta=f_theta(ind) 
     217             dtheta=f_dtheta(ind) 
     218             CALL compute_divgrad(theta,dtheta,1,1) 
    399219          ENDDO 
    400         ENDDO 
    401       ENDDO 
    402  
    403       IF (using_mpi) THEN  
    404         CALL reduce_max_omp(dthetamax ,dthetamax1) 
    405 !$OMP MASTER         
    406         CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    407 !$OMP END MASTER      
    408         CALL bcast_omp(dthetamax)   
    409       ELSE 
    410         CALL allreduce_max_omp(dthetamax,dthetamax1) 
    411         dumax=dumax1 
    412       ENDIF   
    413        
    414       IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 
    415  
    416       DO ind=1,ndomain 
    417         IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    418         CALL swap_dimensions(ind) 
    419         CALL swap_geometry(ind) 
    420         theta=f_theta(ind) 
    421         dtheta=f_dtheta(ind) 
    422         theta=dtheta/dthetamax 
    423       ENDDO 
    424     ENDDO  
    425  
     220       ENDDO        
     221       CALL transfert_request(f_dtheta,req_i1) 
     222       dthetamax = max_scalar(f_dtheta)       
     223       IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 
     224       CALL rescale_field(1./dthetamax, f_dtheta, f_theta) 
     225    END DO 
     226     
     227    cdivgrad=dthetamax**(-1./niterdivgrad)  
    426228    IF (is_master) PRINT *,"divgrad : divgrad",dthetamax 
    427    
    428     cdivgrad=dthetamax**(-1./niterdivgrad)  
    429229    IF (is_master) PRINT *,"cdivgrad : ",cdivgrad 
    430230 
    431       
     231  END SUBROUTINE dissip_constants 
     232 
     233  SUBROUTINE dissip_profile 
     234    USE disvert_mod 
     235    INTEGER    :: l 
     236    REAL(rstd) :: mintau, zz, zvert, fact 
    432237    fact=2 
    433238    DO l=1,llm 
    … …  
    441246       tau_gradrot(l) = tau_gradrot(l)/zvert 
    442247       tau_divgrad(l) = tau_divgrad(l)/zvert 
    443     ENDDO 
    444  
     248    END DO 
     249     
    445250    mintau=tau_graddiv(1) 
    446251    DO l=1,llm 
    447       mintau=MIN(mintau,tau_graddiv(l)) 
    448       mintau=MIN(mintau,tau_gradrot(l)) 
    449       mintau=MIN(mintau,tau_divgrad(l)) 
    450     ENDDO 
    451  
     252       mintau=MIN(mintau,tau_graddiv(l)) 
     253       mintau=MIN(mintau,tau_gradrot(l)) 
     254       mintau=MIN(mintau,tau_divgrad(l)) 
     255    END DO 
     256     
    452257    IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau) 
    453  
     258     
    454259    IF(mintau>0) THEN 
    455260       IF (itau_dissip==0) THEN 
    456          IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..." 
    457          itau_dissip=INT(mintau/dt) 
     261          IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..." 
     262          itau_dissip=INT(mintau/dt) 
    458263       ENDIF 
    459264    ELSE 
    … …  
    464269    dtdissip=itau_dissip*dt 
    465270    IF (is_master) THEN 
    466       PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau 
    467       PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip 
     271       PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau 
     272       PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip 
    468273    ENDIF 
    469  
    470   END SUBROUTINE init_dissip  
    471   
     274     
     275  END SUBROUTINE dissip_profile 
    472276   
    473277  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due) 
    474   USE icosa 
    475   USE theta2theta_rhodz_mod 
    476   USE pression_mod 
    477   USE exner_mod 
    478   USE geopotential_mod 
    479   USE trace 
    480   USE time_mod 
    481   USE omp_para 
    482   IMPLICIT NONE 
    483278    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:) 
    484279    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:) 
    … …  
    492287    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:) 
    493288 
    494     INTEGER :: ind, shear 
    495     INTEGER :: l,ij,nn 
     289    INTEGER :: ind, l,ij,nn 
    496290 
    497291!$OMP BARRIER 
    … …  
    500294    CALL gradiv(f_ue,f_due_diss1) 
    501295    CALL gradrot(f_ue,f_due_diss2) 
    502  
    503296    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss) 
    504297     
    … …  
    516309!DIR$ SIMD 
    517310        DO ij=ij_begin,ij_end 
    518  
    519311            due(ij+u_right,l) = -0.5*( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))*itau_dissip  
    520312            due(ij+u_lup,l)   = -0.5*( due_diss1(ij+u_lup,l)  /tau_graddiv(l) + due_diss2(ij+u_lup,l)  /tau_gradrot(l))*itau_dissip 
    521313            due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip 
    522  
    523314            dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip 
    524315        ENDDO 
    525316      ENDDO 
    526  
    527 !      dtheta_rhodz=0 
    528 !      due=0 
    529  
     317       
    530318      IF(rayleigh_friction_type>0) THEN 
    531        IF(rayleigh_friction_type<99) THEN 
    532          phi=f_geopot(ind) 
    533          ue=f_ue(ind) 
    534          DO l=ll_begin,ll_end 
     319         IF(rayleigh_friction_type<99) THEN 
     320            phi=f_geopot(ind) 
     321            ue=f_ue(ind) 
     322            DO l=ll_begin,ll_end 
     323               DO ij=ij_begin,ij_end 
     324                  CALL relax(t_right, u_right) 
     325                  CALL relax(t_lup,   u_lup) 
     326                  CALL relax(t_ldown, u_ldown) 
     327               ENDDO 
     328            END DO 
     329         ELSE 
     330            ue=f_ue(ind) 
    535331            DO ij=ij_begin,ij_end 
    536                CALL relax(t_right, u_right) 
    537                CALL relax(t_lup,   u_lup) 
    538                CALL relax(t_ldown, u_ldown) 
     332               nn = ij+u_right 
     333               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 
     334                  !print*, "latitude", lat_e(nn)*180./3.14159 
     335                  due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 
     336               ENDIF 
     337               nn = ij+u_lup 
     338               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 
     339                  due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 
     340               ENDIF 
     341               nn = ij+u_ldown 
     342               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 
     343                  due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 
     344               ENDIF 
    539345            ENDDO 
    540          END DO 
    541        ELSE 
    542          ue=f_ue(ind) 
    543             DO ij=ij_begin,ij_end 
    544               nn = ij+u_right 
    545               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 
    546                 !print*, "latitude", lat_e(nn)*180./3.14159 
    547                 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 
    548               ENDIF 
    549               nn = ij+u_lup 
    550               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 
    551                 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 
    552               ENDIF 
    553               nn = ij+u_ldown 
    554               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 
    555                 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 
    556               ENDIF 
    557             ENDDO 
    558        ENDIF 
     346         ENDIF 
    559347      END IF 
    560348   END DO 
    561  
     349    
    562350   CALL trace_end("dissip") 
    563  
     351    
    564352   CALL write_dissip_tendencies 
    565353!$OMP BARRIER 
    566354 
    567     CONTAINS 
    568  
    569       SUBROUTINE relax(shift_t, shift_u) 
    570         USE dcmip_initial_conditions_test_1_2_3 
    571         REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain 
    572              p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain 
    573              fz, u3d(3), uref 
    574         REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values 
    575         LOGICAL :: hybrid_eta 
    576         INTEGER :: shift_u, shift_t, zcoords, nn 
    577         z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g) 
    578         IF(z>zh) THEN  ! relax only in the sponge layer z>zh 
    579  
    580            nn = ij+shift_u 
    581            zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm 
    582            CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, & 
    583                 hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q) 
    584 !           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:) 
    585            u3d = ulon*elon_e(nn,:) ! ulat=0 
    586            uref = sum(u3d*ep_e(nn,:)) 
    587  
    588            fz = sin((pi/2)*(z-zh)/(ztop-zh)) 
    589            fz = fz*fz/rayleigh_tau 
    590            due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref) 
    591 !           fz = 1./rayleigh_tau 
    592 !           due(nn,l) = due(nn,l) - fz*ue(nn,l) 
    593          END IF 
    594        END SUBROUTINE relax 
    595        
    596        SUBROUTINE write_dissip_tendencies 
    597          USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat 
    598          USE wind_mod 
    599          USE output_field_mod 
    600  
    601          CALL transfert_request(f_due_diss1,req_e1_vect) 
    602          CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 
    603          CALL output_field("dulon_diss1",f_buf_ulon) 
    604          CALL output_field("dulat_diss1",f_buf_ulat) 
    605 ! 
    606          CALL transfert_request(f_due_diss2,req_e1_vect) 
    607          CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 
    608          CALL output_field("dulon_diss2",f_buf_ulon) 
    609          CALL output_field("dulat_diss2",f_buf_ulat) 
    610        END SUBROUTINE write_dissip_tendencies 
    611  
    612   END SUBROUTINE dissip 
    613  
    614   SUBROUTINE gradiv(f_ue,f_due) 
    615   USE icosa 
    616   USE trace 
    617   USE omp_para 
    618   IMPLICIT NONE 
    619     TYPE(t_field),POINTER :: f_ue(:) 
    620     TYPE(t_field),POINTER :: f_due(:) 
    621     REAL(rstd),POINTER    :: ue(:,:) 
    622     REAL(rstd),POINTER    :: due(:,:) 
    623     INTEGER :: ind 
    624     INTEGER :: it,l,ij 
    625         
    626     CALL trace_start("gradiv") 
    627  
    628     DO ind=1,ndomain 
     355 CONTAINS 
     356    
     357   SUBROUTINE relax(shift_t, shift_u) 
     358     USE dcmip_initial_conditions_test_1_2_3 
     359     REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain 
     360          p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain 
     361          fz, u3d(3), uref 
     362     REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values 
     363     LOGICAL :: hybrid_eta 
     364     INTEGER :: shift_u, shift_t, zcoords, nn 
     365     z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g) 
     366     IF(z>zh) THEN  ! relax only in the sponge layer z>zh 
     367        nn = ij+shift_u 
     368        zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm 
     369        CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, & 
     370             hyam,hybm,rayleigh_shear,ulon,ulat,w,t,phis,ps,rho,q) 
     371        u3d = ulon*elon_e(nn,:) ! ulat=0 
     372        uref = sum(u3d*ep_e(nn,:)) 
     373         
     374        fz = sin((pi/2)*(z-zh)/(ztop-zh)) 
     375        fz = fz*fz/rayleigh_tau 
     376        due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref) 
     377     END IF 
     378   END SUBROUTINE relax 
     379    
     380   SUBROUTINE write_dissip_tendencies 
     381     USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat 
     382     USE wind_mod 
     383     USE output_field_mod 
     384      
     385     CALL transfert_request(f_due_diss1,req_e1_vect) 
     386     CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 
     387     CALL output_field("dulon_diss1",f_buf_ulon) 
     388     CALL output_field("dulat_diss1",f_buf_ulat) 
     389     ! 
     390     CALL transfert_request(f_due_diss2,req_e1_vect) 
     391     CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 
     392     CALL output_field("dulon_diss2",f_buf_ulon) 
     393     CALL output_field("dulat_diss2",f_buf_ulat) 
     394   END SUBROUTINE write_dissip_tendencies 
     395    
     396 END SUBROUTINE dissip 
     397 
     398  
     399 SUBROUTINE gradiv(f_ue,f_due) 
     400   TYPE(t_field), POINTER :: f_ue(:) 
     401   TYPE(t_field), POINTER :: f_due(:) 
     402   REAL(rstd),POINTER :: ue(:,:) 
     403   REAL(rstd),POINTER :: due(:,:) 
     404   INTEGER :: ind 
     405   INTEGER :: it,l,ij 
     406    
     407   CALL trace_start("gradiv") 
     408    
     409   DO ind=1,ndomain 
    629410      IF (.NOT. assigned_domain(ind)) CYCLE 
    630411      CALL swap_dimensions(ind) 
    … …  
    633414      due=f_due(ind)  
    634415      DO  l = ll_begin, ll_end 
    635 !DIR$ SIMD 
    636         DO ij=ij_begin,ij_end 
    637              due(ij+u_right,l)=ue(ij+u_right,l) 
    638              due(ij+u_lup,l)=ue(ij+u_lup,l) 
    639              due(ij+u_ldown,l)=ue(ij+u_ldown,l) 
    640         ENDDO 
    641       ENDDO 
    642     ENDDO 
    643  
    644     DO it=1,nitergdiv 
    645          
     416         !DIR$ SIMD 
     417         DO ij=ij_begin,ij_end 
     418            due(ij+u_right,l)=ue(ij+u_right,l) 
     419            due(ij+u_lup,l)=ue(ij+u_lup,l) 
     420            due(ij+u_ldown,l)=ue(ij+u_ldown,l) 
     421         ENDDO 
     422      ENDDO 
     423   ENDDO 
     424    
     425   DO it=1,nitergdiv      
    646426      CALL send_message(f_due,req_due) 
    647427      CALL wait_message(req_due) 
    648          
    649428      DO ind=1,ndomain 
    650         IF (.NOT. assigned_domain(ind)) CYCLE 
    651         CALL swap_dimensions(ind) 
    652         CALL swap_geometry(ind) 
    653         due=f_due(ind)  
    654         CALL compute_gradiv(due,due,ll_begin,ll_end) 
    655       ENDDO 
    656     ENDDO 
    657  
     429         IF (.NOT. assigned_domain(ind)) CYCLE 
     430         CALL swap_dimensions(ind) 
     431         CALL swap_geometry(ind) 
     432         due=f_due(ind)  
     433         CALL compute_gradiv(due,due,ll_begin,ll_end) 
     434      ENDDO 
     435   ENDDO 
     436    
    658437   CALL trace_end("gradiv") 
    659  
    660   END SUBROUTINE gradiv 
    661    
    662  
    663   SUBROUTINE gradrot(f_ue,f_due) 
    664   USE icosa 
    665   USE trace 
    666   USE omp_para 
    667   IMPLICIT NONE 
    668     TYPE(t_field),POINTER :: f_ue(:) 
    669     TYPE(t_field),POINTER :: f_due(:) 
    670     REAL(rstd),POINTER    :: ue(:,:) 
    671     REAL(rstd),POINTER    :: due(:,:) 
    672     INTEGER :: ind 
    673     INTEGER :: it,l,ij 
    674         
     438 END SUBROUTINE gradiv 
     439  
     440   
     441 SUBROUTINE gradrot(f_ue,f_due) 
     442   TYPE(t_field), POINTER :: f_ue(:) 
     443    TYPE(t_field), POINTER :: f_due(:) 
     444    REAL(rstd),POINTER :: ue(:,:) 
     445    REAL(rstd),POINTER :: due(:,:) 
     446    INTEGER :: ind, it,l,ij 
    675447    CALL trace_start("gradrot") 
    676448 
    … …  
    692464 
    693465    DO it=1,nitergrot 
    694          
    695       CALL send_message(f_due,req_due) 
    696       CALL wait_message(req_due) 
    697          
    698       DO ind=1,ndomain 
    699         IF (.NOT. assigned_domain(ind)) CYCLE 
    700         CALL swap_dimensions(ind) 
    701         CALL swap_geometry(ind) 
    702         due=f_due(ind)  
    703         CALL compute_gradrot(due,due,ll_begin,ll_end) 
    704       ENDDO 
    705  
    706     ENDDO 
    707  
     466       CALL send_message(f_due,req_due) 
     467       CALL wait_message(req_due) 
     468       DO ind=1,ndomain 
     469          IF (.NOT. assigned_domain(ind)) CYCLE 
     470          CALL swap_dimensions(ind) 
     471          CALL swap_geometry(ind) 
     472          due=f_due(ind)  
     473          CALL compute_gradrot(due,due,ll_begin,ll_end) 
     474       ENDDO 
     475    ENDDO 
     476     
    708477    CALL trace_end("gradrot") 
    709  
    710478  END SUBROUTINE gradrot 
    711479   
    712480  SUBROUTINE divgrad(f_theta,f_dtheta) 
    713   USE icosa 
    714   USE trace 
    715   USE omp_para 
    716   IMPLICIT NONE 
    717     TYPE(t_field),POINTER :: f_theta(:) 
    718     TYPE(t_field),POINTER :: f_dtheta(:) 
    719     REAL(rstd),POINTER    :: theta(:,:) 
    720     REAL(rstd),POINTER    :: dtheta(:,:) 
    721     INTEGER :: ind 
    722     INTEGER :: it 
     481    TYPE(t_field), POINTER :: f_theta(:) 
     482    TYPE(t_field), POINTER :: f_dtheta(:) 
     483    REAL(rstd),POINTER :: theta(:,:) 
     484    REAL(rstd),POINTER :: dtheta(:,:) 
     485    INTEGER :: ind, it 
    723486 
    724487    CALL trace_start("divgrad") 
    … …  
    733496    ENDDO 
    734497 
    735     DO it=1,niterdivgrad 
    736          
    737       CALL transfert_request(f_dtheta,req_i1) 
    738          
     498    DO it=1,niterdivgrad         
     499      CALL transfert_request(f_dtheta,req_i1)         
    739500      DO ind=1,ndomain 
    740501        IF (.NOT. assigned_domain(ind)) CYCLE 
    … …  
    744505        CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end) 
    745506      ENDDO 
    746  
    747507    ENDDO 
    748508 
    749509    CALL trace_end("divgrad") 
    750  
    751510  END SUBROUTINE divgrad 
    752511    
    753512  SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz) 
    754   USE icosa 
    755   USE trace 
    756   USE omp_para 
    757   IMPLICIT NONE 
    758     TYPE(t_field),POINTER :: f_mass(:) 
    759     TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    760     TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
     513    TYPE(t_field), POINTER :: f_mass(:) 
     514    TYPE(t_field), POINTER :: f_theta_rhodz(:) 
     515    TYPE(t_field), POINTER :: f_dtheta_rhodz(:) 
    761516     
    762517    REAL(rstd),POINTER :: mass(:,:) 
    … …  
    813568    ENDDO 
    814569 
    815  
    816570    CALL trace_end("divgrad") 
    817  
    818571  END SUBROUTINE divgrad_theta_rhodz   
    819572   
    820573  SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle) 
    821   USE icosa 
    822   IMPLICIT NONE 
    823574    INTEGER,INTENT(IN)     :: llb 
    824575    INTEGER,INTENT(IN)     :: lle 
    … …  
    843594    DO l=llb,lle 
    844595!DIR$ SIMD 
    845       DO ij=ij_begin,ij_end 
    846   
     596      DO ij=ij_begin,ij_end  
    847597          gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) )        
    848  
    849           gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))        
    850      
     598          gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))     
    851599          gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) ) 
    852  
    853600      ENDDO 
    854601    ENDDO 
    … …  
    867614   
    868615  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle) 
    869   USE icosa 
    870   IMPLICIT NONE 
    871616    INTEGER,INTENT(IN)     :: llb 
    872617    INTEGER,INTENT(IN)     :: lle 
    … …  
    874619    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm) 
    875620    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle) 
    876  
    877     INTEGER :: ij,l 
    878  
    879         
     621    INTEGER :: ij,l        
    880622    DO l=llb,lle 
    881623!DIR$ SIMD 
    882624      DO ij=ij_begin_ext,ij_end_ext 
    883   
    884           grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )         
    885  
    886           grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))        
    887      
     625          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) ) 
     626          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l )) 
    888627          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) ) 
    889  
    890       ENDDO 
    891     ENDDO 
    892      
     628      ENDDO 
     629    ENDDO     
    893630     
    894631    DO l=llb,lle 
    895632!DIR$ SIMD 
    896633      DO ij=ij_begin,ij_end 
    897  
    898634          divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  & 
    899635                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &   
    … …  
    912648 
    913649  END SUBROUTINE compute_divgrad 
    914  
    915650     
    916651  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle) 
    917   USE icosa 
    918   IMPLICIT NONE 
    919652    INTEGER,INTENT(IN)     :: llb 
    920653    INTEGER,INTENT(IN)     :: lle 
    … …  
    943676!DIR$ SIMD 
    944677      DO ij=ij_begin,ij_end 
    945    
    946           gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l))  
    947           
    948           gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l))  
    949          
    950           gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l)) 
    951          
     678          gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l))          
     679          gradrot_e(ij+u_lup,l)  =1/le(ij+u_lup)  *ne(ij,lup)  *(rot_v(ij+z_up,l)   -rot_v(ij+z_lup,l)) 
     680          gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))         
    952681      ENDDO 
    953682    ENDDO 
    … …  
    955684    DO l=llb,lle 
    956685!DIR$ SIMD 
    957       DO ij=ij_begin,ij_end 
    958      
     686      DO ij=ij_begin,ij_end     
    959687          gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot        
    960688          gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot 
    961689          gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot 
    962          
    963690      ENDDO 
    964691    ENDDO   
    … …  
    966693  END SUBROUTINE compute_gradrot 
    967694 
     695!----------------------- Utility routines ------------------ 
     696 
     697  SUBROUTINE global_max(dumax) 
     698    USE mpi_mod 
     699    USE mpipara 
     700    REAL(rstd) :: dumax, dumax1 
     701    IF (using_mpi) THEN  
     702       CALL reduce_max_omp(dumax,dumax1) 
     703       !$OMP MASTER         
     704       CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
     705       !$OMP END MASTER      
     706       CALL bcast_omp(dumax)   
     707    ELSE 
     708       CALL allreduce_max_omp(dumax,dumax1) 
     709       dumax=dumax1 
     710    ENDIF 
     711  END SUBROUTINE global_max     
     712 
     713  FUNCTION max_vel(f_du) 
     714    TYPE(t_field)       :: f_du(:) 
     715    REAL(rstd),POINTER  :: du(:) 
     716    INTEGER :: ind, i,j,ij 
     717    REAL(rstd) :: max_vel, dumax 
     718 
     719    dumax=0. 
     720    DO ind=1,ndomain 
     721       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     722       CALL swap_dimensions(ind) 
     723       CALL swap_geometry(ind) 
     724       du=f_du(ind) 
     725        
     726       DO j=jj_begin,jj_end 
     727          DO i=ii_begin,ii_end 
     728             ij=(j-1)*iim+i    
     729             if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 
     730             if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 
     731             if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 
     732          ENDDO 
     733       ENDDO 
     734    ENDDO 
     735    CALL global_max(dumax) 
     736    max_vel=dumax 
     737  END FUNCTION max_vel 
     738 
     739  FUNCTION max_scalar(f_dtheta) 
     740    TYPE(t_field)       :: f_dtheta(:) 
     741    REAL(rstd),POINTER  :: dtheta(:) 
     742    INTEGER :: ind, i,j,ij 
     743    REAL(rstd) :: max_scalar, dthetamax 
     744    DO ind=1,ndomain 
     745       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     746       CALL swap_dimensions(ind) 
     747       dtheta=f_dtheta(ind) 
     748       DO j=jj_begin,jj_end 
     749          DO i=ii_begin,ii_end 
     750             ij=(j-1)*iim+i 
     751             dthetamax=MAX(dthetamax,ABS(dtheta(ij))) 
     752          ENDDO 
     753       ENDDO 
     754    ENDDO 
     755    CALL global_max(dthetamax) 
     756    max_scalar=dthetamax 
     757  END FUNCTION max_scalar 
     758 
     759  SUBROUTINE random_vel(f_u) 
     760    TYPE(t_field)      :: f_u(:) 
     761    REAL(rstd),POINTER :: u(:) 
     762    INTEGER :: M, ind, i,j,ij 
     763    REAL(rstd) :: r 
     764 
     765    !$OMP BARRIER 
     766    !$OMP MASTER 
     767    DO ind=1,ndomain 
     768       CALL swap_dimensions(ind) 
     769       CALL swap_geometry(ind) 
     770       u=f_u(ind) 
     771        
     772       ! set random seed to get reproductibility when using a different number of process 
     773       CALL RANDOM_SEED(size=M) 
     774       CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
     775        
     776       DO j=jj_begin,jj_end 
     777          DO i=ii_begin,ii_end 
     778             ij=(j-1)*iim+i    
     779             CALL RANDOM_NUMBER(r) 
     780             u(ij+u_right)=r-0.5 
     781             CALL RANDOM_NUMBER(r) 
     782             u(ij+u_lup)=r-0.5 
     783             CALL RANDOM_NUMBER(r) 
     784             u(ij+u_ldown)=r-0.5 
     785          ENDDO 
     786       ENDDO 
     787    ENDDO 
     788    !$OMP END MASTER   
     789    !$OMP BARRIER 
     790     
     791  END SUBROUTINE random_vel 
     792   
     793  SUBROUTINE random_scalar(f_theta) 
     794    TYPE(t_field)      :: f_theta(:) 
     795    REAL(rstd),POINTER :: theta(:) 
     796    INTEGER :: M, ind, i,j,ij 
     797    REAL(rstd) :: r 
     798    !$OMP BARRIER 
     799    !$OMP MASTER 
     800    DO ind=1,ndomain 
     801       CALL swap_dimensions(ind) 
     802       CALL swap_geometry(ind) 
     803       theta=f_theta(ind) 
     804       ! set random seed to get reproductibility when using a different number of process 
     805       CALL RANDOM_SEED(size=M) 
     806       CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
     807        
     808       DO j=jj_begin,jj_end 
     809          DO i=ii_begin,ii_end 
     810             ij=(j-1)*iim+i    
     811             CALL RANDOM_NUMBER(r) 
     812             theta(ij)=r-0.5 
     813          ENDDO 
     814       ENDDO 
     815    ENDDO 
     816    !$OMP END MASTER 
     817    !$OMP BARRIER 
     818  END SUBROUTINE random_scalar 
     819   
     820  SUBROUTINE rescale_field(scal, f_du, f_u) 
     821    TYPE(t_field)      :: f_du(:), f_u(:) 
     822    REAL(rstd),POINTER :: du(:), u(:) 
     823    REAL(rstd) :: scal 
     824    INTEGER :: ind 
     825    DO ind=1,ndomain 
     826       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     827       CALL swap_dimensions(ind) 
     828       u=f_u(ind) 
     829       du=f_du(ind) 
     830       u=scal*du 
     831    ENDDO 
     832  END SUBROUTINE rescale_field 
    968833 
    969834END MODULE dissip_gcm_mod 
    970         
Note: See TracChangeset for help on using the changeset viewer.