da_set_map_para.inc

1 subroutine da_set_map_para
2 
3    !-----------------------------------------------------------------------
4    ! Purpose: TBD
5    !-----------------------------------------------------------------------
6 
7    implicit none
8 
9    real            :: phictr, r, conv
10    real            :: cell, cell2, psx
11 
12    conv = 180.0 / pi
13                                                                   
14    ! define psi1:
15    if (map_projection.eq.1 .or. map_projection.eq.2) then
16       if (phic.lt.0) then 
17          psi1 = 90.+truelat1_3dv
18          psi1 = -psi1
19       else
20         psi1 = 90.-truelat1_3dv
21       end if          
22    else
23       psi1 = 0.
24    end if 
25 
26    psi1 = psi1/conv
27          
28    ! calculate r
29    if (map_projection.ne.3) then
30       psx = (pole-phic)/conv
31       if (map_projection.eq.1) then
32          cell  = earth_radius*sin(psi1)/cone_factor  
33          cell2 = (tan(psx/2.))/(tan(psi1/2.))
34       end if 
35       if (map_projection.eq.2) then
36          cell  = earth_radius*sin(psx)/cone_factor   
37          cell2 = (1. + cos(psi1))/(1. + cos(psx)) 
38       end if  
39       r = cell*(cell2)**cone_factor 
40       ycntr = -r
41    end if  
42 
43    ! for mercator projection, the projection is true at lat at phi1
44    if (map_projection.eq.3) then
45       c2     = earth_radius*cos(psi1) 
46       phictr = phic/conv 
47       cell   = cos(phictr)/(1.0+sin(phictr)) 
48       ycntr  = - c2*log(cell)  
49    end if 
50 
51 end subroutine da_set_map_para
52 
53