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    if (trace_use_dull) call da_trace_entry("da_set_map_para")
13 
14    conv = 180.0 / pi
15                                                                   
16    ! define psi1:
17    if (map_projection==1 .or. map_projection==2) then
18       if (phic.lt.0) then 
19          psi1 = 90.0+truelat1_3dv
20          psi1 = -psi1
21       else
22         psi1 = 90.0-truelat1_3dv
23       end if          
24    else
25       psi1 = 0.0
26    end if 
27 
28    psi1 = psi1/conv
29          
30    ! calculate r
31    if (map_projection.ne.3) then
32       psx = (pole-phic)/conv
33       if (map_projection==1) then
34          cell  = earth_radius*sin(psi1)/cone_factor  
35          cell2 = (tan(psx/2.0))/(tan(psi1/2.0))
36       end if 
37       if (map_projection==2) then
38          cell  = earth_radius*sin(psx)/cone_factor   
39          cell2 = (1.0 + cos(psi1))/(1.0 + cos(psx)) 
40       end if  
41       r = cell*(cell2)**cone_factor 
42       ycntr = -r
43    end if  
44 
45    ! for mercator projection, the projection is true at lat at phi1
46    if (map_projection==3) then
47       c2     = earth_radius*cos(psi1) 
48       phictr = phic/conv 
49       cell   = cos(phictr)/(1.0+sin(phictr)) 
50       ycntr  = - c2*log(cell)  
51    end if 
52 
53    if (trace_use_dull) call da_trace_exit("da_set_map_para")
54 
55 end subroutine da_set_map_para
56 
57