source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/create_limit_unstruct.f90 @ 3895

Last change on this file since 3895 was 3895, checked in by ymipsl, 9 years ago

Make LMDZ5 be compliant to generate initiale state and compute in openMP mode suing unstructured grid.

YM
File size: 10.4 KB
Line 
1MODULE create_limit_unstruct_mod
2    PRIVATE
3    INTEGER, PARAMETER                             :: lmdep=12
4
5    PUBLIC create_limit_unstruct
6
7CONTAINS
8
9  SUBROUTINE create_limit_unstruct
10  USE dimphy
11  USE xios
12  USE ioipsl,             ONLY : ioget_year_len
13  USE time_phylmdz_mod, ONLY : annee_ref
14  USE indice_sol_mod
15  USE phys_state_var_mod
16  USE mod_phys_lmdz_para
17  IMPLICIT NONE
18    INCLUDE "iniprint.h"
19    REAL,    DIMENSION(klon,lmdep)                 :: sic
20    REAL,    DIMENSION(klon,lmdep)                 :: sst
21    REAL,    DIMENSION(klon,lmdep)                 :: rugos
22    REAL,    DIMENSION(klon,lmdep)                 :: albedo
23    REAL,    DIMENSION(klon_mpi,lmdep)             :: sic_mpi
24    REAL,    DIMENSION(klon_mpi,lmdep)             :: sst_mpi
25    REAL,    DIMENSION(klon_mpi,lmdep)             :: rugos_mpi
26    REAL,    DIMENSION(klon_mpi,lmdep)             :: albedo_mpi
27    INTEGER                                        :: ndays
28    REAL                                           :: fi_ice(klon)
29    REAL, ALLOCATABLE                              :: sic_year(:,:)
30    REAL, ALLOCATABLE                              :: sst_year(:,:)
31    REAL, ALLOCATABLE                              :: rugos_year(:,:)
32    REAL, ALLOCATABLE                              :: albedo_year(:,:)
33    REAL, ALLOCATABLE                              :: pctsrf_t(:,:,:)
34    REAL, ALLOCATABLE                              :: phy_bil(:,:)
35    REAL, ALLOCATABLE                              :: sst_year_mpi(:,:)
36    REAL, ALLOCATABLE                              :: rugos_year_mpi(:,:)
37    REAL, ALLOCATABLE                              :: albedo_year_mpi(:,:)
38    REAL, ALLOCATABLE                              :: pctsrf_t_mpi(:,:,:)
39    REAL, ALLOCATABLE                              :: phy_bil_mpi(:,:)
40    INTEGER :: l,k
41    INTEGER :: nbad
42   
43    ndays=ioget_year_len(annee_ref)
44   
45    IF (is_omp_master) THEN
46      CALL xios_recv_field("sic_limit",sic_mpi)
47      CALL xios_recv_field("sst_limit",sst_mpi)
48      CALL xios_recv_field("rugos_limit",rugos_mpi)
49      CALL xios_recv_field("albedo_limit",albedo_mpi)
50    ENDIF
51    CALL scatter_omp(sic_mpi,sic)
52    CALL scatter_omp(sst_mpi,sst)
53    CALL scatter_omp(rugos_mpi,rugos)
54    CALL scatter_omp(albedo_mpi,albedo)
55   
56    ALLOCATE(sic_year(klon,ndays))
57    ALLOCATE(sst_year(klon,ndays))
58    ALLOCATE(rugos_year(klon,ndays))
59    ALLOCATE(albedo_year(klon,ndays))
60    ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
61    ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
62
63! sic
64    CALL time_interpolation(ndays,sic,'gregorian',sic_year)
65    sic_year(:,:)=sic_year(:,:)/100.  ! convert percent to fraction
66    WHERE(sic_year(:,:)>1.0) sic_year(:,:)=1.0    ! Some fractions have some time large negative values
67    WHERE(sic_year(:,:)<0.0) sic_year(:,:)=0.0    ! probably better to apply alse this filter before horizontal interpolation
68   
69! sst
70    CALL time_interpolation(ndays,sst,'gregorian',sst_year)
71    WHERE(sst_year(:,:)<271.38) sst_year(:,:)=271.38
72
73! rugos   
74    DO l=1, lmdep
75      WHERE(NINT(zmasq(:))/=1) rugos(:,l)=0.001
76    ENDDO
77    CALL time_interpolation(ndays,rugos,'360d',rugos_year)
78
79! albedo   
80    CALL time_interpolation(ndays,albedo,'360d',albedo_year)
81
82
83    DO k=1,ndays
84      fi_ice=sic_year(:,k)
85      WHERE(fi_ice>=1.0  ) fi_ice=1.0
86      WHERE(fi_ice<EPSFRA) fi_ice=0.0
87      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
88      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
89
90!!     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
91!!        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
92!!     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
93!!        pctsrf_t(:,is_sic,k)=fi_ice(:)
94!!     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
95        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
96!     END IF
97      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
98      WHERE(1.0-zmasq<EPSFRA)
99        pctsrf_t(:,is_sic,k)=0.0
100        pctsrf_t(:,is_oce,k)=0.0
101      ELSEWHERE
102        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
103          pctsrf_t(:,is_sic,k)=1.0-zmasq
104          pctsrf_t(:,is_oce,k)=0.0
105        ELSEWHERE
106          pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
107          WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
108             pctsrf_t(:,is_oce,k)=0.0
109             pctsrf_t(:,is_sic,k)=1.0-zmasq
110          END WHERE
111        END WHERE
112      END WHERE
113      nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
114      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
115      nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
116      IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
117    END DO
118   
119    ALLOCATE(sst_year_mpi(klon_mpi,ndays))
120    ALLOCATE(rugos_year_mpi(klon_mpi,ndays))
121    ALLOCATE(albedo_year_mpi(klon_mpi,ndays))
122    ALLOCATE(pctsrf_t_mpi(klon_mpi,nbsrf,ndays))
123    ALLOCATE(phy_bil_mpi(klon_mpi,ndays))
124   
125    CALL gather_omp(pctsrf_t   , pctsrf_t_mpi)
126    CALL gather_omp(sst_year   , sst_year_mpi)
127    CALL gather_omp(phy_bil    , phy_bil_mpi)
128    CALL gather_omp(albedo_year, albedo_year_mpi)
129    CALL gather_omp(rugos_year , rugos_year_mpi)
130
131    IF (is_omp_master) THEN
132      CALL xios_send_field("foce_limout",pctsrf_t_mpi(:,is_oce,:))
133      CALL xios_send_field("fsic_limout",pctsrf_t_mpi(:,is_sic,:))
134      CALL xios_send_field("fter_limout",pctsrf_t_mpi(:,is_ter,:))
135      CALL xios_send_field("flic_limout",pctsrf_t_mpi(:,is_lic,:))
136      CALL xios_send_field("sst_limout", sst_year_mpi)
137      CALL xios_send_field("bils_limout",phy_bil_mpi)
138      CALL xios_send_field("alb_limout", albedo_year_mpi)
139      CALL xios_send_field("rug_limout", rugos_year_mpi)
140    ENDIF
141     
142  END SUBROUTINE create_limit_unstruct
143 
144 
145  SUBROUTINE time_interpolation(ndays,field_in,calendar,field_out)
146  USE pchsp_95_m, only: pchsp_95
147  USE pchfe_95_m, only: pchfe_95
148  USE arth_m, only: arth
149  USE dimphy, ONLY : klon
150  USE ioipsl,             ONLY : ioget_year_len
151  USE time_phylmdz_mod, ONLY : annee_ref
152  IMPLICIT NONE
153   INCLUDE "iniprint.h"
154
155   INTEGER,         INTENT(IN)  :: ndays
156   REAL,            INTENT(IN)  :: field_in(klon,lmdep)
157   CHARACTER(LEN=*),INTENT(IN)  :: calendar
158   REAL,            INTENT(OUT) :: field_out(klon,ndays)
159 
160   INTEGER :: ndays_in
161   REAL    :: timeyear(lmdep)   
162   REAL    :: yder(lmdep)   
163   INTEGER :: ij,ierr, n_extrap
164   LOGICAL :: skip
165   
166   ndays_in=year_len(annee_ref, calendar)
167   timeyear=mid_months(annee_ref, calendar, lmdep)
168   
169   n_extrap = 0
170   skip=.FALSE.
171   DO ij=1,klon
172     yder = pchsp_95(timeyear, field_in(ij, :), ibeg=2, iend=2, vc_beg=0., vc_end=0.)
173     CALL pchfe_95(timeyear, field_in(ij, :), yder, skip, arth(0., real(ndays_in) / ndays, ndays), field_out(ij, :), ierr)
174     if (ierr < 0) stop 1
175     n_extrap = n_extrap + ierr
176   END DO
177   
178   IF (n_extrap /= 0) then
179     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
180   ENDIF
181 
182 
183  END SUBROUTINE time_interpolation
184  !-------------------------------------------------------------------------------
185  !
186  FUNCTION year_len(y,cal_in)
187  !
188  !-------------------------------------------------------------------------------
189    USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
190    IMPLICIT NONE
191  !-------------------------------------------------------------------------------
192  ! Arguments:
193    INTEGER                       :: year_len
194    INTEGER,           INTENT(IN) :: y
195    CHARACTER(LEN=*),  INTENT(IN) :: cal_in
196  !-------------------------------------------------------------------------------
197  ! Local variables:
198    CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
199  !-------------------------------------------------------------------------------
200  !--- Getting the input calendar to reset at the end of the function
201    CALL ioget_calendar(cal_out)
202 
203  !--- Unlocking calendar and setting it to wanted one
204    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
205 
206  !--- Getting the number of days in this year
207    year_len=ioget_year_len(y)
208 
209  !--- Back to original calendar
210    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
211 
212  END FUNCTION year_len
213  !
214  !-------------------------------------------------------------------------------
215 
216 
217  !-------------------------------------------------------------------------------
218  !
219  FUNCTION mid_months(y,cal_in,nm)
220  !
221  !-------------------------------------------------------------------------------
222    USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
223    IMPLICIT NONE
224  !-------------------------------------------------------------------------------
225  ! Arguments:
226    INTEGER,                INTENT(IN) :: y               ! year
227    CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
228    INTEGER,                INTENT(IN) :: nm              ! months/year number
229    REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
230  !-------------------------------------------------------------------------------
231  ! Local variables:
232    CHARACTER(LEN=99)                  :: mess            ! error message
233    CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
234    INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
235    INTEGER                            :: m               ! months counter
236    INTEGER                            :: nd              ! number of days
237    INTEGER                            :: k
238  !-------------------------------------------------------------------------------
239    nd=year_len(y,cal_in)
240 
241    IF(nm==12) THEN
242 
243    !--- Getting the input calendar to reset at the end of the function
244      CALL ioget_calendar(cal_out)
245 
246    !--- Unlocking calendar and setting it to wanted one
247      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
248 
249    !--- Getting the length of each month
250      DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
251 
252    !--- Back to original calendar
253      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
254 
255    ELSE IF(MODULO(nd,nm)/=0) THEN
256      WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
257        nm,' months/year. Months number should divide days number.'
258      CALL abort_physic('mid_months',TRIM(mess),1)
259 
260    ELSE
261      mnth=(/(m,m=1,nm,nd/nm)/)
262    END IF
263 
264  !--- Mid-months times
265    mid_months(1)=0.5*REAL(mnth(1))
266    DO k=2,nm
267      mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
268    END DO
269 
270  END FUNCTION mid_months
271 
272
273END MODULE create_limit_unstruct_mod
Note: See TracBrowser for help on using the repository browser.