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

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