Ignore:
Timestamp:
Oct 21, 2024, 2:58:45 PM (23 hours ago)
Author:
abarral
Message:

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File:
1 moved

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/dyn3dmem/friction_loc.F90

    r5245 r5246  
    22! $Id: friction_p.F 1299 2010-01-20 14:27:21Z fairhead $
    33!
    4 c=======================================================================
    5       SUBROUTINE friction_loc(ucov,vcov,pdt)
    6       USE parallel_lmdz
    7       USE control_mod
     4!=======================================================================
     5SUBROUTINE friction_loc(ucov,vcov,pdt)
     6  USE parallel_lmdz
     7  USE control_mod
    88#ifdef CPP_IOIPSL
    9       USE IOIPSL
     9  USE IOIPSL
    1010#else
    11 ! if not using IOIPSL, we still need to use (a local version of) getin
    12       USE ioipsl_getincom
     11  ! if not using IOIPSL, we still need to use (a local version of) getin
     12  USE ioipsl_getincom
    1313#endif
    14       USE comconst_mod, ONLY: pi
    15       IMPLICIT NONE
    16 
    17 !=======================================================================
    18 !
    19 !   Friction for the Newtonian case:
    20 !   --------------------------------
    21 !    2 possibilities (depending on flag 'friction_type'
    22 !    friction_type=0 : A friction that is only applied to the lowermost
    23 !                       atmospheric layer
    24 !    friction_type=1 : Friction applied on all atmospheric layer (but
    25 !       (default)       with stronger magnitude near the surface; see
    26 !                       iniacademic.F)
    27 !=======================================================================
    28 
    29       include "dimensions.h"
    30       include "paramet.h"
    31       include "comgeom2.h"
    32       include "iniprint.h"
    33       include "academic.h"
    34 
    35 ! arguments:
    36       REAL,INTENT(inout) :: ucov( iip1,jjb_u:jje_u,llm )
    37       REAL,INTENT(inout) :: vcov( iip1,jjb_v:jje_v,llm )
    38       REAL,INTENT(in) :: pdt ! time step
    39 
    40 ! local variables:
    41 
    42       REAL modv(iip1,jjb_u:jje_u),zco,zsi
    43       REAL vpn,vps,upoln,upols,vpols,vpoln
    44       REAL u2(iip1,jjb_u:jje_u),v2(iip1,jjb_v:jje_v)
    45       INTEGER i,j,l
    46       REAL,PARAMETER :: cfric=1.e-5
    47       LOGICAL,SAVE :: firstcall=.true.
    48       INTEGER,SAVE :: friction_type=1
    49       CHARACTER(len=20) :: modname="friction_p"
    50       CHARACTER(len=80) :: abort_message
     14  USE comconst_mod, ONLY: pi
     15  IMPLICIT NONE
     16
     17  !=======================================================================
     18  !
     19  !   Friction for the Newtonian case:
     20  !   --------------------------------
     21  !    2 possibilities (depending on flag 'friction_type'
     22  ! friction_type=0 : A friction that is only applied to the lowermost
     23  !                   atmospheric layer
     24  ! friction_type=1 : Friction applied on all atmospheric layer (but
     25  !   (default)       with stronger magnitude near the surface; see
     26  !                   iniacademic.F)
     27  !=======================================================================
     28
     29  include "dimensions.h"
     30  include "paramet.h"
     31  include "comgeom2.h"
     32  include "iniprint.h"
     33  include "academic.h"
     34
     35  ! arguments:
     36  REAL,INTENT(inout) :: ucov( iip1,jjb_u:jje_u,llm )
     37  REAL,INTENT(inout) :: vcov( iip1,jjb_v:jje_v,llm )
     38  REAL,INTENT(in) :: pdt ! time step
     39
     40  ! local variables:
     41
     42  REAL :: modv(iip1,jjb_u:jje_u),zco,zsi
     43  REAL :: vpn,vps,upoln,upols,vpols,vpoln
     44  REAL :: u2(iip1,jjb_u:jje_u),v2(iip1,jjb_v:jje_v)
     45  INTEGER :: i,j,l
     46  REAL,PARAMETER :: cfric=1.e-5
     47  LOGICAL,SAVE :: firstcall=.true.
     48  INTEGER,SAVE :: friction_type=1
     49  CHARACTER(len=20) :: modname="friction_p"
     50  CHARACTER(len=80) :: abort_message
    5151!$OMP THREADPRIVATE(firstcall,friction_type)
    52       integer :: jjb,jje
     52  integer :: jjb,jje
    5353
    5454!$OMP SINGLE
    55       IF (firstcall) THEN
    56         ! set friction type
    57         call getin("friction_type",friction_type)
    58         if ((friction_type.lt.0).or.(friction_type.gt.1)) then
    59           abort_message="wrong friction type"
    60           write(lunout,*)'Friction: wrong friction type',friction_type
    61           call abort_gcm(modname,abort_message,42)
    62         endif
    63         firstcall=.false.
    64       ENDIF
     55  IF (firstcall) THEN
     56    ! ! set friction type
     57    call getin("friction_type",friction_type)
     58    if ((friction_type.lt.0).or.(friction_type.gt.1)) then
     59      abort_message="wrong friction type"
     60      write(lunout,*)'Friction: wrong friction type',friction_type
     61      call abort_gcm(modname,abort_message,42)
     62    endif
     63    firstcall=.false.
     64  ENDIF
    6565!$OMP END SINGLE COPYPRIVATE(friction_type,firstcall)
    6666
    67       if (friction_type.eq.0) then ! friction on first layer only
     67  if (friction_type.eq.0) then ! friction on first layer only
    6868!$OMP SINGLE
    69 c   calcul des composantes au carre du vent naturel
    70       jjb=jj_begin
    71       jje=jj_end+1
    72       if (pole_sud) jje=jj_end
    73      
    74       do j=jjb,jje
    75          do i=1,iip1
    76             u2(i,j)=ucov(i,j,1)*ucov(i,j,1)*unscu2(i,j)
    77          enddo
    78       enddo
    79      
    80       jjb=jj_begin-1
    81       jje=jj_end+1
    82       if (pole_nord) jjb=jj_begin
    83       if (pole_sud) jje=jj_end-1
    84      
    85       do j=jjb,jje
    86          do i=1,iip1
    87             v2(i,j)=vcov(i,j,1)*vcov(i,j,1)*unscv2(i,j)
    88          enddo
    89       enddo
    90 
    91 c   calcul du module de V en dehors des poles
    92       jjb=jj_begin
    93       jje=jj_end+1
    94       if (pole_nord) jjb=jj_begin+1
    95       if (pole_sud) jje=jj_end-1
    96      
    97       do j=jjb,jje
    98          do i=2,iip1
    99             modv(i,j)=sqrt(0.5*(u2(i-1,j)+u2(i,j)+v2(i,j-1)+v2(i,j)))
    100          enddo
    101          modv(1,j)=modv(iip1,j)
    102       enddo
    103 
    104 c   les deux composantes du vent au pole sont obtenues comme
    105 c   premiers modes de fourier de v pres du pole
    106       if (pole_nord) then
    107      
    108         upoln=0.
    109         vpoln=0.
    110      
    111         do i=2,iip1
    112            zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
    113            zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
    114            vpn=vcov(i,1,1)/cv(i,1)
    115            upoln=upoln+zco*vpn
    116            vpoln=vpoln+zsi*vpn
    117         enddo
    118         vpn=sqrt(upoln*upoln+vpoln*vpoln)/pi
    119         do i=1,iip1
    120 c          modv(i,1)=vpn
    121            modv(i,1)=modv(i,2)
    122         enddo
    123 
    124       endif
    125      
    126       if (pole_sud) then
    127      
    128         upols=0.
    129         vpols=0.
    130         do i=2,iip1
    131            zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
    132            zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
    133            vps=vcov(i,jjm,1)/cv(i,jjm)
    134            upols=upols+zco*vps
    135            vpols=vpols+zsi*vps
    136         enddo
    137         vps=sqrt(upols*upols+vpols*vpols)/pi
    138         do i=1,iip1
    139 c        modv(i,jjp1)=vps
    140          modv(i,jjp1)=modv(i,jjm)
    141         enddo
    142      
    143       endif
    144      
    145 c   calcul du frottement au sol.
    146 
    147       jjb=jj_begin
    148       jje=jj_end
    149       if (pole_nord) jjb=jj_begin+1
    150       if (pole_sud) jje=jj_end-1
    151 
    152       do j=jjb,jje
    153          do i=1,iim
    154             ucov(i,j,1)=ucov(i,j,1)
    155      s      -cfric*pdt*0.5*(modv(i+1,j)+modv(i,j))*ucov(i,j,1)
    156          enddo
    157          ucov(iip1,j,1)=ucov(1,j,1)
    158       enddo
    159      
    160       jjb=jj_begin
    161       jje=jj_end
    162       if (pole_sud) jje=jj_end-1
    163      
    164       do j=jjb,jje
    165          do i=1,iip1
    166             vcov(i,j,1)=vcov(i,j,1)
    167      s      -cfric*pdt*0.5*(modv(i,j+1)+modv(i,j))*vcov(i,j,1)
    168          enddo
    169          vcov(iip1,j,1)=vcov(1,j,1)
    170       enddo
     69  !   calcul des composantes au carre du vent naturel
     70  jjb=jj_begin
     71  jje=jj_end+1
     72  if (pole_sud) jje=jj_end
     73
     74  do j=jjb,jje
     75     do i=1,iip1
     76        u2(i,j)=ucov(i,j,1)*ucov(i,j,1)*unscu2(i,j)
     77     enddo
     78  enddo
     79
     80  jjb=jj_begin-1
     81  jje=jj_end+1
     82  if (pole_nord) jjb=jj_begin
     83  if (pole_sud) jje=jj_end-1
     84
     85  do j=jjb,jje
     86     do i=1,iip1
     87        v2(i,j)=vcov(i,j,1)*vcov(i,j,1)*unscv2(i,j)
     88     enddo
     89  enddo
     90
     91  !   calcul du module de V en dehors des poles
     92  jjb=jj_begin
     93  jje=jj_end+1
     94  if (pole_nord) jjb=jj_begin+1
     95  if (pole_sud) jje=jj_end-1
     96
     97  do j=jjb,jje
     98     do i=2,iip1
     99        modv(i,j)=sqrt(0.5*(u2(i-1,j)+u2(i,j)+v2(i,j-1)+v2(i,j)))
     100     enddo
     101     modv(1,j)=modv(iip1,j)
     102  enddo
     103
     104  !   les deux composantes du vent au pole sont obtenues comme
     105  !   premiers modes de fourier de v pres du pole
     106  if (pole_nord) then
     107
     108    upoln=0.
     109    vpoln=0.
     110
     111    do i=2,iip1
     112       zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
     113       zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
     114       vpn=vcov(i,1,1)/cv(i,1)
     115       upoln=upoln+zco*vpn
     116       vpoln=vpoln+zsi*vpn
     117    enddo
     118    vpn=sqrt(upoln*upoln+vpoln*vpoln)/pi
     119    do i=1,iip1
     120       ! modv(i,1)=vpn
     121       modv(i,1)=modv(i,2)
     122    enddo
     123
     124  endif
     125
     126  if (pole_sud) then
     127
     128    upols=0.
     129    vpols=0.
     130    do i=2,iip1
     131       zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
     132       zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
     133       vps=vcov(i,jjm,1)/cv(i,jjm)
     134       upols=upols+zco*vps
     135       vpols=vpols+zsi*vps
     136    enddo
     137    vps=sqrt(upols*upols+vpols*vpols)/pi
     138    do i=1,iip1
     139     ! modv(i,jjp1)=vps
     140     modv(i,jjp1)=modv(i,jjm)
     141    enddo
     142
     143  endif
     144
     145  !   calcul du frottement au sol.
     146
     147  jjb=jj_begin
     148  jje=jj_end
     149  if (pole_nord) jjb=jj_begin+1
     150  if (pole_sud) jje=jj_end-1
     151
     152  do j=jjb,jje
     153     do i=1,iim
     154        ucov(i,j,1)=ucov(i,j,1) &
     155              -cfric*pdt*0.5*(modv(i+1,j)+modv(i,j))*ucov(i,j,1)
     156     enddo
     157     ucov(iip1,j,1)=ucov(1,j,1)
     158  enddo
     159
     160  jjb=jj_begin
     161  jje=jj_end
     162  if (pole_sud) jje=jj_end-1
     163
     164  do j=jjb,jje
     165     do i=1,iip1
     166        vcov(i,j,1)=vcov(i,j,1) &
     167              -cfric*pdt*0.5*(modv(i,j+1)+modv(i,j))*vcov(i,j,1)
     168     enddo
     169     vcov(iip1,j,1)=vcov(1,j,1)
     170  enddo
    171171!$OMP END SINGLE
    172       endif ! of if (friction_type.eq.0)
    173 
    174       if (friction_type.eq.1) then
    175        ! for ucov()
    176         jjb=jj_begin
    177         jje=jj_end
    178         if (pole_nord) jjb=jj_begin+1
    179         if (pole_sud) jje=jj_end-1
     172  endif ! of if (friction_type.eq.0)
     173
     174  if (friction_type.eq.1) then
     175   ! ! for ucov()
     176    jjb=jj_begin
     177    jje=jj_end
     178    if (pole_nord) jjb=jj_begin+1
     179    if (pole_sud) jje=jj_end-1
    180180
    181181!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    182         do l=1,llm
    183           ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)*
    184      &                                  (1.-pdt*kfrict(l))
    185         enddo
     182    do l=1,llm
     183      ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)* &
     184            (1.-pdt*kfrict(l))
     185    enddo
    186186!$OMP END DO NOWAIT
    187        
    188        ! for vcoc()
    189         jjb=jj_begin
    190         jje=jj_end
    191         if (pole_sud) jje=jj_end-1
    192        
     187
     188   ! ! for vcoc()
     189    jjb=jj_begin
     190    jje=jj_end
     191    if (pole_sud) jje=jj_end-1
     192
    193193!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    194         do l=1,llm
    195           vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)*
    196      &                                  (1.-pdt*kfrict(l))
    197         enddo
     194    do l=1,llm
     195      vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)* &
     196            (1.-pdt*kfrict(l))
     197    enddo
    198198!$OMP END DO
    199       endif ! of if (friction_type.eq.1)
    200 
    201       RETURN
    202       END
    203 
     199  endif ! of if (friction_type.eq.1)
     200
     201  RETURN
     202END SUBROUTINE friction_loc
     203
Note: See TracChangeset for help on using the changeset viewer.