source: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_limit_unstruct.F90 @ 5437

Last change on this file since 5437 was 3323, checked in by Laurent Fairhead, 7 years ago

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